stdErrLiar = function(morVal, datafr){
  prediction <- datafr %>%
    filter(role == "Liar", mor == morVal) %>%
    mutate(probabilityRed = p,
           drawnRed = ks,
           predicted = val) %>%
    select(expt, probabilityRed, drawnRed, predicted)
  real <- bs.final %>%
    filter(roleCurrent == "bullshitter") %>%
    select(expt, probabilityRed, drawnRed, reportedDrawn)
  left_join(real, prediction, by=c("expt", "probabilityRed", "drawnRed")) %>%
    mutate(error = reportedDrawn - predicted) %>%
    summarise(n = n(),
              stderror = sum(error^2)/sqrt(n)) %>%
    .$stderror
}

stdErrDetector = function(morVal, datafr){
  prediction <- datafr %>%
    filter(role == "Detector", mor == morVal) %>%
    mutate(probabilityRed = p,
           reportedDrawn = ks,
           predicted = val) %>%
    select(expt, probabilityRed, reportedDrawn, predicted)
  real <- bs.final %>%
    filter(roleCurrent == "bullshitDetector") %>%
    select(expt, probabilityRed, reportedDrawn, callBS) %>%
    mutate(callBS = ifelse(callBS=="True", 1, 0))
  left_join(real, prediction, by=c("expt", "probabilityRed", "reportedDrawn")) %>%
    mutate(error = callBS - predicted) %>%
    summarise(n = n(),
              stderror = sum(error^2)/sqrt(n)) %>%
    .$stderror
}

data.frame(mor = min(df.unif$mor):max(df.unif$mor)) %>%
  mutate(stderrLieMorality = sapply(mor, stdErrLiar, df.unif),
         stderrDetectMorality = sapply(mor, stdErrDetector, df.unif),
         stderrLieMoralityIdiot = sapply(mor, stdErrLiar, moralityResponseToIdiot)) %>%
  drop_na() %>%
  arrange(stderrLieMoralityIdiot) %>%
  filter(stderrLieMorality == min(stderrLieMorality) | 
           stderrDetectMorality == min(stderrDetectMorality) | 
           stderrLieMoralityIdiot == min(stderrLieMoralityIdiot))
##   mor stderrLieMorality stderrDetectMorality stderrLieMoralityIdiot
## 1  18          230.8268             22.07332               233.5609
## 2  20          230.4797             22.50151               236.3627
## 3   8          312.9083             20.53429               598.3536
lieMoral = 20
detectMoral = 8

Pretty Plots

Recursive Inference Model

# Recursive Model w/ Morality Penalty

lie.recurse.moral <- df.unif %>%
  filter(role=="Liar", mor==lieMoral, expt=="expt4") %>%
  mutate(p = paste("p =", p)) %>%
  ggplot(aes(x=ks, y=val, colour=p)) +
  geom_abline(size=ablineSize, intercept = 0, slope = 1, colour="darkgray", linetype=2) +
  geom_line(size=lineSize, alpha=alphaLine) +
  ggtitle("Recursive+Aversion") +
  scale_x_continuous("", labels=rep("",6), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
lie.recurse.moral

ggsave("img/modelRecurse_lie_moral.png", lie.recurse.moral, height=2.5, width=3)

detect.recurse.moral <- df.unif %>%
  filter(role=="Detector", mor==detectMoral, expt=="expt4") %>%
  mutate(p = paste("p =", p)) %>%
  ggplot(aes(x=ks, y=val, colour=p)) +
  geom_line(size=lineSize, alpha=alphaLine) +
  scale_x_continuous("", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
detect.recurse.moral

ggsave("img/modelRecurse_detect_moral.png", detect.recurse.moral, height=8, width=7.5)
lie.recurse.nomoral <- df.unif %>%
  filter(role=="Liar", mor==0, expt=="expt4") %>%
  mutate(p = paste("p =", p)) %>%
  ggplot(aes(x=ks, y=val, colour=p)) +
  geom_abline(size=ablineSize, intercept = 0, slope = 1, colour="darkgray", linetype=2) +
  geom_line(size=lineSize, alpha=alphaLine) +
  ggtitle("Recursive") +
  scale_x_continuous("", labels=rep("",6), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
lie.recurse.nomoral

ggsave("img/modelRecurse_lie_nomoral.png", lie.recurse.nomoral, height=8, width=7.5)

detect.recurse.nomoral <- df.unif %>%
  filter(role=="Detector", mor==0, expt=="expt4") %>%
  mutate(p = paste("p =", p)) %>%
  ggplot(aes(x=ks, y=val, colour=p)) +
  geom_line(size=lineSize, alpha=alphaLine) +
  scale_x_continuous("", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
detect.recurse.nomoral

ggsave("img/modelRecurse_detect_nomoral.png", detect.recurse.nomoral, height=8, width=7.5)

Response to Idiot Liar/Lie Detector

moral = lieMoral
Expt = 4
LIE <- 1-diag(numMarbles+1)
mat.idiotLiar <- matrix(rep(1/11,11*11), nrow=11)
P.K <- matrix(rep(p.k(0:10, 0.2), each=10+1), nrow=10+1)
lie0.2idiot <- rowSums(P.K*mat.idiotLiar*LIE)/rowSums(P.K*mat.idiotLiar)
detect0.2idiot <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.2, lie0.2idiot)

P.K <- matrix(rep(p.k(0:10, 0.5), each=10+1), nrow=10+1)
lie0.5idiot <- rowSums(P.K*mat.idiotLiar*LIE)/rowSums(P.K*mat.idiotLiar)
detect0.5idiot <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.5, lie0.5idiot)

P.K <- matrix(rep(p.k(0:10, 0.8), each=10+1), nrow=10+1)
lie0.8idiot <- rowSums(P.K*mat.idiotLiar*LIE)/rowSums(P.K*mat.idiotLiar)
detect0.8idiot <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.8, lie0.8idiot)

detect.respondToIdiot <- data.frame(p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
           ksay=rep(0:10,3),
           bs=c(detect0.2idiot, detect0.5idiot, detect0.8idiot)) %>%
  mutate(probabilityRed.txt = paste("p =", p)) %>%
  ggplot(aes(x=ksay, y=bs, colour=probabilityRed.txt)) +
  geom_line(size=lineSize, alpha=alphaLine) +
  # scale_x_continuous("", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  # scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_x_continuous("Reported Marbles", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("Proportion BS Called", limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff)+
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.text.y = element_text(margin=margin(0,-1,0,-2.5)),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
detect.respondToIdiot

#ggsave("img/modelRespondToIdiot_detect.png", detect.respondToIdiot, height=8, width=7.5)           

avgPropBS <- bs.final %>%
  filter(roleCurrent == "bullshitDetector") %>%
  summarise(propBS = sum(callBS == "True")/n()) %>%
  .$propBS
bonddepaulo <- 1-.56
vrij <- 1 - .615
lie.respondToIdiot <- data.frame(p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
           k=rep(0:10,3),
           ksay=c(exp.ksay(0.2, rep(bonddepaulo, 11)),
                  exp.ksay(0.5, rep(bonddepaulo, 11)),
                  exp.ksay(0.8, rep(bonddepaulo, 11)))) %>%
  mutate(probabilityRed.txt = paste("p =", p)) %>%
  ggplot(aes(x=k, y=ksay, colour=probabilityRed.txt)) +
  geom_abline(intercept = 0, slope = 1, colour="darkgray", linetype=2, size=ablineSize) +
  geom_line(size=lineSize, alpha=alphaLine) +
  # ggtitle("Uniform + Moral") +
  ggtitle("Uniform Opponent") +
  # scale_x_continuous("", labels=rep("",6), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  # scale_y_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_x_continuous("True Marbles", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("Reported Marbles", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.text.y = element_text(margin=margin(0,-1,0,3.8)),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
lie.respondToIdiot

#ggsave("img/modelRespondToIdiot_lie.png", lie.respondToIdiot, height=8, width=7.5)  

Response to Idiot Liar/Lie Detector Without Morality

moral = 0
Expt = 4
detect0.2idiot.nomoral <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.2, lie0.2idiot)
detect0.5idiot.nomoral <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.5, lie0.5idiot)
detect0.8idiot.nomoral <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.8, lie0.8idiot)

detect.respondToIdiot.nomoral <- data.frame(p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
           ksay=rep(0:10,3),
           bs=c(detect0.2idiot, detect0.5idiot, detect0.8idiot)) %>%
  mutate(probabilityRed.txt = paste("p =", p)) %>%
  ggplot(aes(x=ksay, y=bs, colour=probabilityRed.txt)) +
  geom_line(size=lineSize, alpha=alphaLine) +
  # scale_x_continuous("", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  # scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_x_continuous("Reported Marbles", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("Proportion BS Called", limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff)+
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.text.y = element_text(margin=margin(0,-1,0,-2.5)),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
detect.respondToIdiot.nomoral

ggsave("img/modelRespondToIdiot_detect_nomoral.png", detect.respondToIdiot.nomoral, height=8, width=7.5)           

lie.respondToIdiot.nomoral <- data.frame(p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
           k=rep(0:10,3),
           ksay=c(exp.ksay(0.2, rep(bonddepaulo, 11)),
                  exp.ksay(0.5, rep(bonddepaulo, 11)),
                  exp.ksay(0.8, rep(bonddepaulo, 11)))) %>%
  mutate(probabilityRed.txt = paste("p =", p)) %>%
  ggplot(aes(x=k, y=ksay, colour=probabilityRed.txt)) +
  geom_abline(intercept = 0, slope = 1, colour="darkgray", linetype=2, size=ablineSize) +
  geom_line(size=lineSize, alpha=alphaLine) +
  # ggtitle("Uniform + Moral") +
  ggtitle("Random Opponent") +
  # scale_x_continuous("", labels=rep("",6), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  # scale_y_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_x_continuous("True Marbles", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("Reported Marbles", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.text.y = element_text(margin=margin(0,-1,0,3.8)),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
lie.respondToIdiot.nomoral

ggsave("img/modelRespondToIdiot_lie_nomoral.png", lie.respondToIdiot.nomoral, height=8, width=7.5)  
moral = lieMoral
lie.results_expt4 <- bs.final %>%
  filter(roleCurrent == "bullshitter", expt=="expt4") %>%
  mutate(probabilityRed.txt = paste("p =", probabilityRed)) %>%
  group_by(probabilityRed.txt, drawnRed) %>%
  summarise(mean = mean(reportedDrawn),
            sd = sd(reportedDrawn),
            se = sd/sqrt(n()),
            n=n()) %>%
  filter(n > 3) %>%
  ggplot(aes(x=drawnRed, y=mean, colour=probabilityRed.txt, fill=probabilityRed.txt)) +
  geom_abline(intercept = 0, slope = 1, colour="darkgray", linetype=2, size=ablineSize) +
  geom_ribbon(aes(ymin=mean-se, ymax=mean+se), alpha=0.3) +
  geom_line(size=0.8, alpha=alphaLine) +
  geom_point(pch=21, stroke=0.5, size=1.5, colour="black", alpha=alphaLine) +
  ggtitle("Red Utility") +
  #scale_x_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  #scale_y_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_colour_manual(values = condColors.diff) +
  scale_fill_manual(values = condColors.diff) +
  guides(fill=guide_legend(title="", override.aes=list(size=2.5)), colour="none") +
  theme_minimal() +
  theme(legend.position="none",
        legend.text = element_text(size=7),
        legend.key.size = unit(0.36, 'cm'),
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
lie.results_expt4

ggsave("img/results_lie_expt4.png", lie.results_expt4, height=2.5, width=3)

detect.results_expt4 <- bs.final %>%
  filter(roleCurrent == "bullshitDetector", expt=="expt4") %>%
  mutate(probabilityRed.txt = paste("p =", probabilityRed)) %>%
  group_by(probabilityRed.txt, expt, reportedDrawn) %>%
  summarise(prop = (sum(callBS == "True")+1)/(n()+2),
            se = sqrt((prop*(1-prop)/n())),
            n = n()) %>%
  filter(n > 3) %>%
  ggplot(aes(x=reportedDrawn, y=prop, colour=probabilityRed.txt, fill=probabilityRed.txt)) +
  geom_ribbon(aes(ymin=prop-se, ymax=prop+se), alpha=0.3) +
  geom_line(size=0.8, alpha=alphaLine) +
  geom_point(pch=21, stroke=0.5, size=1.5, colour="black", alpha=alphaLine) +
  #scale_x_continuous("", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  #scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_x_continuous("Reported Marbles", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("Proportion BS Called", limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values = condColors.diff, guide="none") +
  scale_fill_manual(values = condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
detect.results_expt4

ggsave("img/results_detect_expt4.png", detect.results_expt4, height=2.5, width=3)
lie.results_expt5 <- bs.final %>%
  filter(roleCurrent == "bullshitter", expt=="expt5") %>%
  mutate(probabilityRed.txt = paste("p =", probabilityRed)) %>%
  group_by(probabilityRed.txt, drawnRed) %>%
  summarise(mean = mean(reportedDrawn),
            sd = sd(reportedDrawn),
            se = sd/sqrt(n()),
            n=n()) %>%
  filter(n > 3) %>%
  ggplot(aes(x=drawnRed, y=mean, colour=probabilityRed.txt, fill=probabilityRed.txt)) +
  geom_abline(intercept = 0, slope = 1, colour="darkgray", linetype=2, size=ablineSize) +
  geom_ribbon(aes(ymin=mean-se, ymax=mean+se), alpha=0.3) +
  geom_line(size=0.8, alpha=alphaLine) +
  geom_point(pch=21, stroke=0.5, size=1.5, colour="black", alpha=alphaLine) +
  ggtitle("Blue Utility") +
  scale_x_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_colour_manual(values = condColors.diff) +
  scale_fill_manual(values = condColors.diff) +
  guides(fill=guide_legend(title="", override.aes=list(size=2.5)), colour="none") +
  theme_minimal() +
  theme(legend.text = element_text(size=7),
        legend.key.size = unit(0.36, 'cm'),
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
lie.results_expt5 + theme(legend.position="none")

ggsave("img/results_lie_expt5.png", lie.results_expt5, height=2.5, width=3)

detect.results_expt5 <- bs.final %>%
  filter(roleCurrent == "bullshitDetector", expt=="expt5") %>%
  mutate(probabilityRed.txt = paste("p =", probabilityRed)) %>%
  group_by(probabilityRed.txt, reportedDrawn) %>%
  summarise(prop = (sum(callBS == "True")+1)/(n()+2),
            se = sqrt((prop*(1-prop)/n())),
            n = n()) %>%
  filter(n > 3) %>%
  ggplot(aes(x=reportedDrawn, y=prop, colour=probabilityRed.txt, fill=probabilityRed.txt)) +
  geom_ribbon(aes(ymin=prop-se, ymax=prop+se), alpha=0.3) +
  geom_line(size=0.8, alpha=alphaLine) +
  geom_point(pch=21, stroke=0.5, size=1.5, colour="black", alpha=alphaLine) +
  scale_x_continuous("", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values = condColors.diff, guide="none") +
  scale_fill_manual(values = condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
detect.results_expt5

ggsave("img/results_detect_expt5.png", detect.results_expt5, height=2.5, width=3)

Combine Graphs

legend <- get_legend(lie.results_expt5)

lieplots <- plot_grid(#lie.idiot,
                      lie.respondToIdiot.nomoral,
                      #lie.respondToIdiot,
                      #lie.respondToHumanData,
                      lie.recurse.nomoral,
                      lie.recurse.moral,
                      lie.results_expt4,
                      lie.results_expt5 + theme(legend.position="none"), 
                      nrow=1,
                      rel_widths=c(10.92, 10, 10, 10, 10)) #10.72 with 4 panels

detectplots <- plot_grid(#detect.idiot, 
                         detect.respondToIdiot.nomoral,
                         #detect.respondToIdiot, 
                         #detect.respondToHumanData, 
                         detect.recurse.nomoral, 
                         detect.recurse.moral, 
                         detect.results_expt4,
                         detect.results_expt5, 
                         nrow=1,
                         rel_widths=c(10.92, 10, 10, 10, 10))

liarLabel <- ggdraw() +
  draw_label("Sender", fontface="bold", size=12, y=0.87, x=0.5, hjust=0, angle=270)

lierow <- plot_grid(lieplots, liarLabel, nrow=1, rel_widths=c(96, 4))

detectorLabel <- ggdraw() +
  draw_label("Receiver", fontface="bold", size=12, y=0.99, x=0.5, hjust=0, angle=270)

detectrow <- plot_grid(detectplots, detectorLabel, nrow=1, rel_widths=c(96, 4))

allplots <- plot_grid(lierow, 
                      detectrow, 
                      nrow=2,
                      #rel_heights=c(0.53, 0.47))
                      rel_heights=c(0.538, 0.462))

withLegend <- ggdraw() +
  draw_plot(allplots) +
  draw_plot(legend, x=0.415, y=0.16)

# ggsave("img/allPredictions_wExpt5.pdf", 
#        withLegend, 
#        units="cm", height=8.5, width=17.8)

withLegend

ggsave("img/allPredictions_wExpt5.pdf", 
       withLegend, 
       units="cm", height=7.3, width=17.8)

# ggsave("img/allPredictions_wExpt5_nolegend.png", 
#        allplots, 
#        units="cm", height=7.3, width=17.8)


detectplots_nomoral <- plot_grid(detect.respondToIdiot.nomoral + ggtitle("Random Opponent") + theme(plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold")),
                         detect.recurse.nomoral + ggtitle("Recursive") + theme(plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold")), 
                         detect.results_expt4 + ggtitle("Red Utility") + theme(plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold")),
                         detect.results_expt5 + ggtitle("Blue Utility") + theme(plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold")), 
                         nrow=1,
                         rel_widths=c(10.72, 10, 10, 10))


detectrow_nomoral <- plot_grid(detectplots_nomoral, detectorLabel, nrow=1, rel_widths=c(96, 4))

withLegend_nomoral <- ggdraw() +
  draw_plot(detectrow_nomoral) +
  draw_plot(legend, x=0.415, y=0.32) 
ggsave("img/detectPredictions.pdf", 
       withLegend_nomoral, 
       units="cm", height=4.7, width=17.8)






detectplots_results <- plot_grid( detect.results_expt4 + ggtitle("Red Response") + theme(plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold")),
                         detect.results_expt5 + ggtitle("Blue Response") + theme(plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold")), 
                         nrow=1,
                         rel_widths=c(11.32, 10))

detectorLabel <- ggdraw() +
  draw_label("Receiver", fontface="bold", size=12, y=0.85, x=0.5, hjust=0, angle=270)

detectrow_results <- plot_grid(detectplots_results, detectorLabel, nrow=1, rel_widths=c(94, 6))

withLegend_results <- ggdraw() +
  draw_plot(detectrow_results) +
  draw_plot(legend, x=0.35, y=0.31) 
ggsave("img/detectResults.pdf", 
       withLegend_results, 
       units="cm", height=4.3, width=8.7)
bs.final %>%
  filter(roleCurrent == "bullshitter", expt=="expt4") %>%
  mutate(probabilityRed.txt = paste("p =", probabilityRed)) %>%
  group_by(probabilityRed.txt, drawnRed, reportedDrawn) %>%
  summarise(n=n())
## # A tibble: 246 x 4
## # Groups:   probabilityRed.txt, drawnRed [33]
##    probabilityRed.txt drawnRed reportedDrawn     n
##    <chr>                 <int>         <int> <int>
##  1 p = 0.2                   0             0    60
##  2 p = 0.2                   0             1    15
##  3 p = 0.2                   0             2    26
##  4 p = 0.2                   0             3    21
##  5 p = 0.2                   0             4    10
##  6 p = 0.2                   0             5     6
##  7 p = 0.2                   0             8     1
##  8 p = 0.2                   0             9     1
##  9 p = 0.2                   0            10     3
## 10 p = 0.2                   1             0    14
## # … with 236 more rows
p <- plot_ly(z = ~volcano) %>% add_surface()

Model vs Human Comparison

modelColors = c("random opponent"="chocolate2",
                "uniform+aversion"="firebrick3",
                "recursive "="goldenrod1", 
                "recursive+aversion"="forestgreen") #breaks if no space after "recursive"

pShapes = c("p = 0.2" = 21,
            "p = 0.5" = 22,
            "p = 0.8" = 24
)

prob_to_logit <- function(p){
  log(p/(1-p))
}

modelLieComp <- df.unif %>%
  filter(role=="Liar", mor==lieMoral) %>%
  mutate(model = "recursive+aversion", 
         k = ks,
         mean.model = val,
         se.model = se) %>%
  select(role, model, p, k, expt, mean.model, se.model)

modelDetectComp <- df.unif %>%
  filter(role=="Detector", mor==detectMoral) %>%
  mutate(model = "recursive+aversion", 
         k = ks,
         mean.model = val,
         se.model = se) %>%
  select(role, model, p, k, expt, mean.model, se.model)

rationalLieComp <- df.unif %>%
  filter(role=="Liar", mor==0) %>%
  mutate(model = "recursive ", 
         k = ks,
         mean.model = val,
         se.model = se) %>%
  select(role, model, p, k, expt, mean.model, se.model)

rationalDetectComp <- df.unif %>%
  filter(role=="Detector", mor==0) %>%
  mutate(model = "recursive ", 
         k = ks,
         mean.model = val,
         se.model = se) %>%
  select(role, model, p, k, expt, mean.model, se.model)

Expt = 4
moral = lieMoral
idiotLieComp <- data.frame(role = "Liar",
                           model = "uniform+aversion",
                           p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                           k=rep(0:10,3),
                           expt = "expt4",
                           mean.model=c(exp.ksay(0.2, rep(bonddepaulo, 11)),
                                  exp.ksay(0.5, rep(bonddepaulo, 11)),
                                  exp.ksay(0.8, rep(bonddepaulo, 11))),
                           se.model = NA)

idiotDetectComp <- data.frame(role = "Detector",
                              model = "uniform+aversion",
                              p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                              k=rep(0:10,3),
                              expt = "expt4",
                              mean.model=c(detect0.2idiot, 
                                           detect0.5idiot, 
                                           detect0.8idiot),
                              se.model = NA)

moral = 0
idiotLieComp.nomoral <- data.frame(role = "Liar",
                           model = "random opponent",
                           p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                           k=rep(0:10,3),
                           expt = "expt4",
                           mean.model=c(exp.ksay(0.2, rep(bonddepaulo, 11)),
                                  exp.ksay(0.5, rep(bonddepaulo, 11)),
                                  exp.ksay(0.8, rep(bonddepaulo, 11))),
                           se.model = NA)

idiotDetectComp.nomoral <- data.frame(role = "Detector",
                              model = "random opponent",
                              p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                              k=rep(0:10,3),
                              expt = "expt4",
                              mean.model=c(detect0.2idiot.nomoral, 
                                           detect0.5idiot.nomoral, 
                                           detect0.8idiot.nomoral),
                              se.model = NA)


Expt = 5
moral = lieMoral
P.K <- matrix(rep(p.k(0:10, 0.2), each=10+1), nrow=10+1)
lie0.2idiot <- rowSums(P.K*mat.idiotLiar*LIE)/rowSums(P.K*mat.idiotLiar)
detect0.2idiot_expt5 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.2, lie0.2idiot)

P.K <- matrix(rep(p.k(0:10, 0.5), each=10+1), nrow=10+1)
lie0.5idiot <- rowSums(P.K*mat.idiotLiar*LIE)/rowSums(P.K*mat.idiotLiar)
detect0.5idiot_expt5 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.5, lie0.5idiot)

P.K <- matrix(rep(p.k(0:10, 0.8), each=10+1), nrow=10+1)
lie0.8idiot <- rowSums(P.K*mat.idiotLiar*LIE)/rowSums(P.K*mat.idiotLiar)
detect0.8idiot_expt5 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.8, lie0.8idiot)

idiotLieComp_expt5 <- data.frame(role = "Liar",
                           model = "uniform+aversion",
                           p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                           k=rep(0:10,3),
                           expt = "expt5",
                           mean.model=c(exp.ksay(0.2, rep(bonddepaulo, 11)),
                                  exp.ksay(0.5, rep(bonddepaulo, 11)),
                                  exp.ksay(0.8, rep(bonddepaulo, 11))),
                           se.model = NA)

idiotDetectComp_expt5 <- data.frame(role = "Detector",
                              model = "uniform+aversion",
                              p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                              k=rep(0:10,3),
                              expt = "expt5",
                              mean.model=c(detect0.2idiot_expt5,
                                           detect0.5idiot_expt5,
                                           detect0.8idiot_expt5),
                              se.model = NA)

moral = 0
detect0.2idiot.nomoral_expt5 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.2, lie0.2idiot)
detect0.5idiot.nomoral_expt5 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.5, lie0.5idiot)
detect0.8idiot.nomoral_expt5 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.8, lie0.8idiot)

idiotLieComp.nomoral_expt5 <- data.frame(role = "Liar",
                           model = "random opponent",
                           p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                           k=rep(0:10,3),
                           expt = "expt5",
                           mean.model=c(exp.ksay(0.2, rep(bonddepaulo, 11)),
                                  exp.ksay(0.5, rep(bonddepaulo, 11)),
                                  exp.ksay(0.8, rep(bonddepaulo, 11))),
                           se.model = NA)

idiotDetectComp.nomoral_expt5 <- data.frame(role = "Detector",
                              model = "random opponent",
                              p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                              k=rep(0:10,3),
                              expt = "expt5",
                              mean.model=c(detect0.2idiot.nomoral_expt5,
                                           detect0.5idiot.nomoral_expt5,
                                           detect0.8idiot.nomoral_expt5),
                              se.model = NA)

idiotLieComp <- bind_rows(idiotLieComp, idiotLieComp_expt5)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
idiotDetectComp <- bind_rows(idiotDetectComp, idiotDetectComp_expt5)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
idiotLieComp.nomoral <- bind_rows(idiotLieComp.nomoral, idiotLieComp.nomoral_expt5)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
idiotDetectComp.nomoral <- bind_rows(idiotDetectComp.nomoral, idiotDetectComp.nomoral_expt5)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
modelComp <- bind_rows(idiotLieComp, idiotDetectComp, idiotLieComp.nomoral, idiotDetectComp.nomoral, modelLieComp, modelDetectComp, rationalLieComp, rationalDetectComp)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
humanCompLie <- bs.final %>%
  filter(roleCurrent == "bullshitter") %>%
  mutate(k = drawnRed,
         p = probabilityRed) %>%
  group_by(p, k, expt, roleCurrent) %>%
  summarise(mean.human = mean(reportedDrawn), #not sure if this is kosher, but putting lying on same scale as detecting
            sd = sd(reportedDrawn),
            se.human = sd/sqrt(n()),
            n=n()) %>%
  mutate(role = "Liar") %>%
  select(role, p, k, expt, mean.human, se.human)
humanCompDetect <- bs.final %>%
  filter(roleCurrent == "bullshitDetector") %>%
  mutate(k = reportedDrawn,
         p = probabilityRed) %>%
  group_by(p, k, expt, roleCurrent) %>%
  summarise(mean.human = (sum(callBS == "True")+1)/(n()+2),
            se.human = sqrt((mean.human*(1-mean.human)/n())),
            n=n()) %>%
  mutate(role = "Detector") %>%
  select(role, p, k, expt, mean.human, se.human)
humanComp <- bind_rows(humanCompLie, humanCompDetect)


modelVsHuman <- left_join(modelComp, humanComp, by=c("role", "expt", "p", "k")) %>%
  mutate(model = factor(model, levels=c("random opponent", "uniform+aversion", "recursive ", "recursive+aversion")),
         mean.model.logit = ifelse(role=="Detector", prob_to_logit(mean.model), mean.model),
         mean.human.logit = ifelse(role=="Detector", prob_to_logit(mean.human), mean.human))
## Warning: Column `expt` joining character vector and factor, coercing into
## character vector
## Warning in log(p/(1 - p)): NaNs produced

## Warning in log(p/(1 - p)): NaNs produced
modelVsHuman %>%
  arrange(p, role, k, expt) %>%
  head(10)
##        role              model   p k  expt mean.model   se.model mean.human
## 1  Detector   uniform+aversion 0.2 0 expt4  0.4675204         NA  0.2500000
## 2  Detector    random opponent 0.2 0 expt4  0.4675204         NA  0.2500000
## 3  Detector recursive+aversion 0.2 0 expt4  0.2781886 0.02007997  0.2500000
## 4  Detector         recursive  0.2 0 expt4  0.4036800 0.02407225  0.2500000
## 5  Detector   uniform+aversion 0.2 0 expt5  0.9872513         NA  0.6600791
## 6  Detector    random opponent 0.2 0 expt5  0.9872513         NA  0.6600791
## 7  Detector recursive+aversion 0.2 0 expt5  0.5475090 0.03503614  0.6600791
## 8  Detector         recursive  0.2 0 expt5  0.5945403 0.03409133  0.6600791
## 9  Detector   uniform+aversion 0.2 1 expt4  0.5230766         NA  0.2011494
## 10 Detector    random opponent 0.2 1 expt4  0.5230766         NA  0.2011494
##      se.human mean.model.logit mean.human.logit
## 1  0.05685735      -0.13010146       -1.0986123
## 2  0.05685735      -0.13010146       -1.0986123
## 3  0.05685735      -0.95346471       -1.0986123
## 4  0.05685735      -0.39015494       -1.0986123
## 5  0.02989855       4.34949272        0.6636465
## 6  0.02989855       4.34949272        0.6636465
## 7  0.02989855       0.19061112        0.6636465
## 8  0.02989855       0.38276708        0.6636465
## 9  0.03056525       0.09237218       -1.3791259
## 10 0.03056525       0.09237218       -1.3791259
lieComp <- modelVsHuman %>%
  filter(role=="Liar", model!="uniform+aversion", expt=="expt4") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, fill=model, shape=p.txt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=3) +
  #geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.1) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.1) +
  scale_x_continuous("Model", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0.025,0)) +
  scale_y_continuous("Human", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0.025,0)) +
  scale_fill_manual(values=modelColors) +
  scale_shape_manual(values=pShapes) +
  guides(shape=guide_legend(order=1, title=""),
         fill=guide_legend(order=0, title="", override.aes=list(shape=21))) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=10),
        axis.text = element_text(size=9),
        axis.text.y = element_text(margin=margin(0,0,0,5)),
        axis.line = element_line(size=1, colour=axisLineColour))

detectComp <- modelVsHuman %>%
  filter(role=="Detector", model!="uniform+aversion", expt=="expt4") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, fill=model, shape=p.txt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=3) +
  #geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.02) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.02) +
  scale_x_continuous("Model", limits=c(0,1), expand=c(0,0,0.025,0)) +
  scale_y_continuous("Human", limits=c(0,1), expand=c(0,0,0.025,0)) +
  scale_fill_manual(values=modelColors) +
  scale_shape_manual(values=pShapes) +
  # guides(shape=guide_legend(title="priors"),
  #        fill=guide_legend(override.aes=list(shape=21))) +
  guides(shape=guide_legend(order=1, title=""),
         fill=guide_legend(order=0, title="", override.aes=list(shape=21))) +
  theme_minimal() +
  theme(legend.key.size = unit(0.34, 'cm'),
        legend.text = element_text(size=7),
        legend.box = "horizontal",
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=10),
        axis.text = element_text(size=9),
        axis.text.y = element_text(margin=margin(0,-0.5,0,-2)),
        axis.line = element_line(size=1, colour=axisLineColour))

legendComp <- get_legend(detectComp)
lieCompLabel <- ggdraw() +
  draw_label("Liar", fontface="bold", size=12, y=0.95, x=0.4, hjust=0, angle=270)

detectCompLabel <- ggdraw() + 
  draw_label("Detector", fontface="bold", size=12, y=0.95, x=0.4, hjust=0, angle=270)

bothComp <- plot_grid(lieComp, lieCompLabel, 
                      detectComp + theme(legend.position = "none"), detectCompLabel, ncol=2, rel_widths=c(95,5))

comparison <- ggdraw() +
  draw_plot(bothComp) +
  draw_plot(legendComp, x=-0.08, y=0.45)
  #draw_plot(legendComp, x=-0.18, y=0.35)
comparison

ggsave("img/modelComparison.pdf", comparison, units="cm", height=15, width=8.7)

Model Comparison Facetted by Model

lieCompFacet <- modelVsHuman %>%
  filter(role=="Liar", model!="uniform+aversion", expt=="expt4") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.1) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.1) +
  ggtitle("Lying") +
  scale_x_continuous("Model", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0.025,0)) +
  scale_y_continuous("Human", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0.025,0)) +
  scale_colour_manual(values=condColors.diff) +
  guides(colour=guide_legend(title="", override.aes=list(size=2.5))) +
  facet_wrap(~model) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=10),
        axis.text = element_text(size=9),
        axis.text.y = element_text(margin=margin(0,0,0,5)),
        axis.line = element_line(size=1, colour=axisLineColour))
lieCompFacet

ggsave("img/lieCompareFacet.png", lieCompFacet, height=4, width=13)

detectCompFacet <- modelVsHuman %>%
  filter(role=="Detector", model!="uniform+aversion", expt=="expt4") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.02) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.02) +
  ggtitle("Detecting") +
  scale_x_continuous("Model", limits=c(0,1), expand=c(0,0,0.025,0)) +
  scale_y_continuous("Human", limits=c(0,1), expand=c(0,0,0.025,0)) +
  scale_colour_manual(values=condColors.diff, guide="none") +
  facet_wrap(~model) +
  theme_minimal() +
  theme(legend.key.size = unit(0.34, 'cm'),
        legend.text = element_text(size=7),
        legend.box = "horizontal",
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=10),
        axis.text = element_text(size=9),
        axis.text.y = element_text(margin=margin(0,-0.5,0,-2)),
        axis.line = element_line(size=1, colour=axisLineColour))
detectCompFacet

ggsave("img/detectCompareFacet.png", detectCompFacet, height=4, width=13)





plotTitleSizeComp = 12
axisTitleSizeComp = 10
axisTextSizeComp = 9

shapeCompare <- c(
  "expt4" = 19,
  "expt5" = 1
)

shapeCompare2 <- c(
  "red utility" = 19,
  "blue utility" = 1
)

#Broken down by graphs
lieCompUnif <- modelVsHuman %>%
  filter(role=="Liar", model=="random opponent") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt, shape=expt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=2, alpha=0.9) +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.1) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.1) +
  ggtitle("Random Opponent") +
  scale_x_continuous("Model", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0.025,0)) +
  scale_y_continuous("Human", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0.025,0)) +
  scale_colour_manual(values=condColors.diff) +
  scale_shape_manual(values=shapeCompare) +
  guides(colour=guide_legend(title="", override.aes=list(size=2.5))) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(size=plotTitleSizeComp, face = "bold"),
        plot.margin = marginDim,
        axis.title = element_text(size=axisTitleSizeComp),
        axis.text = element_text(size=axisTextSizeComp),
        axis.text.y = element_text(margin=margin(0,0,0,5)),
        axis.line = element_line(size=1, colour=axisLineColour))
lieCompUnif

#ggsave("img/lieCompareUnif.png", lieCompUnif, height=4.1, width=4)

lieCompRecurse <- modelVsHuman %>%
  filter(role=="Liar", model=="recursive ") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt, shape=expt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=2, alpha=0.9) +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.1) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.1) +
  ggtitle("Recursive") +
  scale_x_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0.025,0)) +
  scale_y_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0.025,0)) +
  scale_colour_manual(values=condColors.diff) +
  scale_shape_manual(values=shapeCompare) +
  guides(colour=guide_legend(title="", override.aes=list(size=2.5))) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(size=plotTitleSizeComp, face = "bold"),
        plot.margin = marginDim,
        axis.title = element_text(size=axisTitleSizeComp),
        axis.text = element_text(size=axisTextSizeComp),
        axis.line = element_line(size=1, colour=axisLineColour))
lieCompRecurse

#ggsave("img/lieCompareRecurse.png", lieCompRecurse, height=4.1, width=4)

lieCompAverse <- modelVsHuman %>%
  filter(role=="Liar", model=="recursive+aversion") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt, shape=expt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=2, alpha=0.9) +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.1) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.1) +
  ggtitle("Recursive + Aversion") +
  scale_x_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0.025,0)) +
  scale_y_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0.025,0)) +
  scale_colour_manual(values=condColors.diff) +
  scale_shape_manual(values=shapeCompare) +
  guides(colour=guide_legend(title="", override.aes=list(size=2.5))) +
  theme_minimal() +
  theme(legend.text = element_text(size=7),
        legend.key.size = unit(0.36, 'cm'),
        plot.title = element_text(size=plotTitleSizeComp, face = "bold"),
        plot.margin = marginDim,
        axis.title = element_text(size=axisTitleSizeComp),
        axis.text = element_text(size=axisTextSizeComp),
        axis.line = element_line(size=1, colour=axisLineColour))
lieCompAverse

#ggsave("img/lieCompareAverse.png", lieCompAverse, height=4.1, width=4)

detectCompUnif <- modelVsHuman %>%
  filter(role=="Detector", model=="random opponent") %>%
  mutate(p.txt = paste("p =", p),
         expt = factor(ifelse(expt == "expt4", "red utility", "blue utility"), levels=c("red utility", "blue utility"))) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt, shape=expt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=2, alpha=0.9) +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.02) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.02) +
  scale_x_continuous("Model", limits=c(0,1), expand=c(0,0,0.025,0)) +
  scale_y_continuous("Human", limits=c(0,1), expand=c(0,0,0.025,0)) +
  scale_colour_manual(values=condColors.diff, guide="none") +
  scale_shape_manual(values=shapeCompare2, guide="none") +
  guides(shape=guide_legend(title="", override.aes=list(size=2.5))) +
  theme_minimal() +
  theme(legend.key.size = unit(0.34, 'cm'),
        legend.text = element_text(size=7),
        legend.text.align = 1,
        legend.box = "horizontal",
        plot.title = element_text(size=plotTitleSize, face = "bold"),
        plot.margin = marginDim,
        axis.title = element_text(size=axisTitleSizeComp),
        axis.text = element_text(size=axisTextSizeComp),
        axis.text.y = element_text(margin=margin(0,-0.5,0,-2)),
        axis.line = element_line(size=1, colour=axisLineColour))
detectCompUnif

#ggsave("img/detectCompareUnif.png", detectCompUnif, height=4.1, width=4)

detectCompRecurse <- modelVsHuman %>%
  filter(role=="Detector", model=="recursive ") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt, shape=expt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=2, alpha=0.9) +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.02) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.02) +
  scale_x_continuous("", limits=c(0,1), expand=c(0,0,0.025,0)) +
  scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0,0.025,0)) +
  scale_colour_manual(values=condColors.diff, guide="none") +
  scale_shape_manual(values=shapeCompare, guide="none") +
  theme_minimal() +
  theme(legend.key.size = unit(0.34, 'cm'),
        legend.text = element_text(size=7),
        legend.box = "horizontal",
        plot.title = element_text(size=plotTitleSize, face = "bold"),
        plot.margin = marginDim,
        axis.title = element_text(size=axisTitleSizeComp),
        axis.text = element_text(size=axisTextSizeComp),
        axis.line = element_line(size=1, colour=axisLineColour))
detectCompRecurse

#ggsave("img/detectCompareRecurse.png", detectCompRecurse, height=4.1, width=4)

detectCompAverse <- modelVsHuman %>%
  filter(role=="Detector", model=="recursive+aversion") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt, shape=expt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=2, alpha=0.9) +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.02) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.02) +
  scale_x_continuous("", limits=c(0,1), expand=c(0,0,0.025,0)) +
  scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0,0.025,0)) +
  scale_colour_manual(values=condColors.diff, guide="none") +
  scale_shape_manual(values=shapeCompare, guide="none") +
  theme_minimal() +
  theme(legend.key.size = unit(0.34, 'cm'),
        legend.text = element_text(size=7),
        legend.box = "horizontal",
        plot.title = element_text(size=plotTitleSize, face = "bold"),
        plot.margin = marginDim,
        axis.title = element_text(size=axisTitleSizeComp),
        axis.text = element_text(size=axisTextSizeComp),
        axis.line = element_line(size=1, colour=axisLineColour))
detectCompAverse

#ggsave("img/detectCompareAverse.png", detectCompAverse, height=4.1, width=4)


legendComp <- get_legend(lie.results_expt5)
legendComp2 <- get_legend(detectCompUnif)

lieCompPlots <- plot_grid(lieCompUnif,
                          lieCompRecurse,
                          lieCompAverse + theme(legend.position="none"), 
                          nrow=1,
                          rel_widths=c(10.9, 10, 10))

detectCompPlots <- plot_grid(detectCompUnif + theme(legend.position="none"),
                             detectCompRecurse,
                             detectCompAverse, 
                             nrow=1,
                             rel_widths=c(10.9, 10, 10))

liarCompLabel <- ggdraw() +
  draw_label("Sender", fontface="bold", size=12, y=0.87, x=0.5, hjust=0, angle=270)

lieComprow <- plot_grid(lieCompPlots, liarCompLabel, nrow=1, rel_widths=c(96, 4))

detectorCompLabel <- ggdraw() +
  draw_label("Receiver", fontface="bold", size=12, y=0.99, x=0.5, hjust=0, angle=270)

detectComprow <- plot_grid(detectCompPlots, detectorCompLabel, nrow=1, rel_widths=c(96, 4))

allCompplots <- plot_grid(lieComprow, 
                      detectComprow, 
                      nrow=2,
                      rel_heights=c(0.525, 0.475))

CompwithLegend <- ggdraw() +
  draw_plot(allCompplots) +
  draw_plot(legendComp, x=0.41, y=0.16) +
  draw_plot(legendComp2, x=0.41, y=0.25) +
  draw_label(expression(paste("r"^"2"*" = 0.118")), x=0.135, y=0.93, size=11) +
  draw_label(expression(paste("r"^"2"*" = 0.477")), x=0.45, y=0.93, size=11) +
  draw_label(expression(paste("r"^"2"*" = 0.965")), x=0.76, y=0.93, size=11) +
  draw_label(expression(paste("r"^"2"*" = 0.237")), x=0.135, y=0.455, size=11) +
  draw_label(expression(paste("r"^"2"*" = 0.596")), x=0.45, y=0.455, size=11) +
  draw_label(expression(paste("r"^"2"*" = 0.730")), x=0.76, y=0.455, size=11)

CompwithLegend

ggsave("img/allComparisons.pdf", 
       CompwithLegend, 
       units="cm", height=12.2, width=17.8)

Compute R^2

# by role and model and expt
modelVsHuman %>%
  group_by(role, model, expt) %>%
  summarise(r = cor(mean.model, mean.human),
            r.squared = r^2,
            p = cor.test(mean.model, mean.human)$p.value)
## # A tibble: 16 x 6
## # Groups:   role, model [8]
##    role     model              expt      r r.squared        p
##    <chr>    <fct>              <chr> <dbl>     <dbl>    <dbl>
##  1 Detector random opponent    expt4 0.501     0.251 2.97e- 3
##  2 Detector random opponent    expt5 0.473     0.224 5.42e- 3
##  3 Detector uniform+aversion   expt4 0.501     0.251 2.97e- 3
##  4 Detector uniform+aversion   expt5 0.473     0.224 5.42e- 3
##  5 Detector "recursive "       expt4 0.758     0.574 3.31e- 7
##  6 Detector "recursive "       expt5 0.796     0.633 3.14e- 8
##  7 Detector recursive+aversion expt4 0.833     0.694 1.76e- 9
##  8 Detector recursive+aversion expt5 0.884     0.781 9.37e-12
##  9 Liar     random opponent    expt4 0.730     0.534 1.39e- 6
## 10 Liar     random opponent    expt5 0.665     0.442 2.46e- 5
## 11 Liar     uniform+aversion   expt4 0.993     0.987 1.24e-30
## 12 Liar     uniform+aversion   expt5 0.976     0.952 5.16e-22
## 13 Liar     "recursive "       expt4 0.886     0.784 7.44e-12
## 14 Liar     "recursive "       expt5 0.848     0.720 4.55e-10
## 15 Liar     recursive+aversion expt4 0.993     0.985 5.05e-30
## 16 Liar     recursive+aversion expt5 0.976     0.953 4.12e-22
# by role and model
modelVsHuman %>%
  group_by(role, model) %>%
  summarise(r = cor(mean.model, mean.human),
            r.squared = r^2,
            p = cor.test(mean.model, mean.human)$p.value)
## # A tibble: 8 x 5
## # Groups:   role [2]
##   role     model                  r r.squared        p
##   <chr>    <fct>              <dbl>     <dbl>    <dbl>
## 1 Detector random opponent    0.487     0.237 3.42e- 5
## 2 Detector uniform+aversion   0.487     0.237 3.42e- 5
## 3 Detector "recursive "       0.772     0.596 3.17e-14
## 4 Detector recursive+aversion 0.854     0.730 7.42e-20
## 5 Liar     random opponent    0.344     0.118 4.65e- 3
## 6 Liar     uniform+aversion   0.982     0.964 6.36e-48
## 7 Liar     "recursive "       0.690     0.477 1.42e-10
## 8 Liar     recursive+aversion 0.982     0.965 3.56e-48
# collapsing over role
modelVsHuman %>%
  group_by(model) %>%
  summarise(r = cor(mean.model, mean.human),
            r.squared = r^2,
            p = cor.test(mean.model, mean.human)$p.value)
## # A tibble: 4 x 4
##   model                  r r.squared         p
##   <fct>              <dbl>     <dbl>     <dbl>
## 1 random opponent    0.741     0.550 2.79e- 24
## 2 uniform+aversion   0.981     0.962 6.32e- 94
## 3 "recursive "       0.896     0.803 9.54e- 48
## 4 recursive+aversion 0.991     0.981 1.65e-114

Model Comparison with Panels for Each Model

lieComp_facet <- modelVsHuman %>%
  filter(model != "uniform+moral", expt == "expt4", role == "Liar") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=1) +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.1) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.1) +
  scale_x_continuous("Model", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("Human", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  facet_grid(model~.) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=10),
        axis.text = element_text(size=9),
        axis.text.y = element_text(margin=margin(0,0,0,5))) +
  annotate("segment", x=-Inf, xend=-Inf, y=-Inf, yend=Inf, colour=axisLineColour, size=2) +
  annotate("segment", x=-Inf, xend=Inf, y=-Inf, yend=-Inf, colour=axisLineColour, size=2)
ggsave("img/lieComparison_facet.png", lieComp_facet, height=7.75, width=3.5)


detectComp_facet <- modelVsHuman %>%
  filter(model != "uniform+moral", expt == "expt4", role == "Detector") %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt)) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  geom_point(size=1) +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.02) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.02) +
  scale_x_continuous("Model", limits=c(0,1), expand=c(0,0)) +
  scale_y_continuous("Human", limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  # guides(colour=guide_legend(override.aes=list(shape=21))) +
  guides(colour=guide_legend(title="")) +
  facet_grid(model~.) +
  theme_minimal() +
  theme(legend.position = "none",
        #legend.key.size = unit(0.34, 'cm'),
        #legend.text = element_text(size=7),
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=10),
        axis.text = element_text(size=9),
        axis.text.y = element_text(margin=margin(0,-0.5,0,-2))) +
  annotate("segment", x=-Inf, xend=-Inf, y=-Inf, yend=Inf, colour=axisLineColour, size=2) +
  annotate("segment", x=-Inf, xend=Inf, y=-Inf, yend=-Inf, colour=axisLineColour, size=2)
ggsave("img/detectComparison_facet.png", detectComp_facet, height=7.75, width=3.5)

Model Predictions Feeding in Human Results

Lie Detector Behavior from Lying Results

moral = 20
Expt = 4

# "template" for tidyr complete() function
template <- data.frame(expt = NA,
                       probabilityRed = rep(0,11),
                       drawnRed = rep(0,11),
                       reportedDrawn = 0:10,
                       n = rep(0,11))

liePropsDF <- bs.final %>%
  filter(roleCurrent == "bullshitter") %>%
  group_by(expt, probabilityRed, drawnRed, reportedDrawn) %>%
  summarise(n=n()) 

liePropsDF <- bind_rows(template,liePropsDF) %>%
  complete(expt, probabilityRed, drawnRed, reportedDrawn, fill=list(n=0)) %>%
  filter(!is.na(expt), probabilityRed != 0) %>%
  group_by(expt, probabilityRed, drawnRed) %>%
  mutate(total = sum(n),
         proportion = n / total) %>%
  filter(expt == "expt4")

LIE <- 1-diag(numMarbles+1)
# equivalent to output of p.L_ksay.k.r
df0.2 <- liePropsDF %>%
  filter(probabilityRed == 0.2)
mat0.2 <- matrix(df0.2$proportion, nrow=11)
# output p_t.ksay.r
P.K <- matrix(rep(p.k(0:numMarbles, 0.2), each=numMarbles+1), nrow=numMarbles+1)
lie0.2 <- rowSums(P.K*mat0.2*LIE)/rowSums(P.K*mat0.2)
# output p.D_bs.ksay.r
detect0.2 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.2, lie0.2)

df0.5 <- liePropsDF %>%
  filter(probabilityRed == 0.5)
mat0.5 <- matrix(df0.5$proportion, nrow=11)
P.K <- matrix(rep(p.k(0:10, 0.5), each=10+1), nrow=10+1)
lie0.5 <- rowSums(P.K*mat0.5*LIE)/rowSums(P.K*mat0.5)
detect0.5 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.5, lie0.5)

df0.8 <- liePropsDF %>%
  filter(probabilityRed == 0.8)
mat0.8 <- matrix(df0.8$proportion, nrow=11)
P.K <- matrix(rep(p.k(0:10, 0.8), each=10+1), nrow=10+1)
lie0.8 <- rowSums(P.K*mat0.8*LIE)/rowSums(P.K*mat0.8)
detect0.8 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.8, lie0.8)

detect.respondToHumanData <- data.frame(p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
           ksay=rep(0:10,3),
           bs=c(detect0.2, detect0.5, detect0.8)) %>%
  mutate(probabilityRed.txt = paste("p =", p)) %>%
  ggplot(aes(x=ksay, y=bs, colour=probabilityRed.txt)) +
  geom_line(size=lineSize, alpha=alphaLine) +
  scale_x_continuous("", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff)+
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
detect.respondToHumanData 

ggsave("img/modelRespondToData_detect.png", detect.respondToHumanData, height=8, width=7.5)

Liar Behavior from Lie Detector Results

bsPropsDF <- bs.final %>%
  filter(roleCurrent == "bullshitDetector", expt=="expt4") %>%
  group_by(probabilityRed, reportedDrawn) %>%
  summarise(count = sum(callBS == "True"),
            total = n(),
            propBS = count/total) %>%
  select(probabilityRed, reportedDrawn, propBS) %>%
  mutate(proptruth = p_t.ksay.r(probabilityRed, propBS),
         lie = exp.ksay(probabilityRed, propBS))
bsProps0.2 <- bsPropsDF %>%
  filter(probabilityRed == 0.2)
propksay0.2 <- colSums(KSAY*p.L_ksay.k.r(0.2, bsProps0.2$propBS))
proptruth <- p_t.ksay.r(0.2, bsProps0.2$propBS)

lie.respondToHumanData <- bsPropsDF %>%
  mutate(probabilityRed.txt = paste("p =", probabilityRed)) %>%
  ggplot(aes(x=reportedDrawn, y=lie, colour=probabilityRed.txt)) +
  geom_abline(intercept = 0, slope = 1, colour="darkgray", linetype=2, size=ablineSize) +
  geom_line(size=lineSize, alpha=alphaLine) +
  ggtitle("Input: Human") +
  scale_x_continuous("", labels=rep("",6), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
lie.respondToHumanData

ggsave("img/modelRespondToData_lie.png", lie.respondToHumanData, height=8, width=7.5)

Uniform Liar + Detector

lie.idiot <- data.frame(p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
           k=rep(0:10,3),
           ksay=sum(0:10 * 1/11)) %>%
  mutate(probabilityRed.txt = paste("p =", p)) %>%
  ggplot(aes(x=k, y=ksay, colour=probabilityRed.txt)) +
  geom_abline(intercept = 0, slope = 1, colour="darkgray", linetype=2, size=ablineSize) +
  geom_line(size=lineSize, alpha=alphaLine) +
  ggtitle("Uniform Opponent") +
  scale_x_continuous("True Marbles", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("Reported Marbles", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.y = element_text(margin=margin(0,0,0,25.5)),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
lie.idiot

ggsave("img/modelIdiot_lie.png", lie.idiot, height=8, width=7.5)     

detect.idiot <- data.frame(p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
           ksay=rep(0:10,3),
           bs=avgPropBS) %>%
  mutate(probabilityRed.txt = paste("p =", p)) %>%
  ggplot(aes(x=ksay, y=bs, colour=probabilityRed.txt)) +
  geom_line(size=lineSize, alpha=alphaLine) +
  scale_x_continuous("Reported Marbles", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("Proportion BS Called", limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff)+
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
detect.idiot

ggsave("img/modelIdiot_detect.png", detect.idiot, height=8, width=7.5)     
# with just morality model
modelLieComp <- df.unif %>%
  filter(role=="Liar", mor==lieMoral) %>%
  mutate(model = "moral", 
         k = ks,
         mean.model = val/10,
         se.model = se/10) %>%
  select(role, model, p, k , mean.model, se.model)

modelDetectComp <- df.unif %>%
  filter(role=="Detector", mor==detectMoral) %>%
  mutate(model = "moral", 
         k = ks,
         mean.model = val,
         se.model = se) %>%
  select(role, model, p, k , mean.model, se.model)

modelComp <- bind_rows(modelLieComp, modelDetectComp)

humanCompLie <- bs.final %>%
  filter(roleCurrent == "bullshitter") %>%
  mutate(k = drawnRed,
         p = probabilityRed) %>%
  group_by(p, k, roleCurrent) %>%
  summarise(mean.human = mean(reportedDrawn)/10, #not sure if this is kosher, but putting lying on same scale as detecting
            sd = sd(reportedDrawn)/10,
            se.human = sd/sqrt(n()),
            n=n()) %>%
  mutate(role = "Liar") %>%
  select(role, p, k , mean.human, se.human)
humanCompDetect <- bs.final %>%
  filter(roleCurrent == "bullshitDetector") %>%
  mutate(k = reportedDrawn,
         p = probabilityRed) %>%
  group_by(p, k, roleCurrent) %>%
  summarise(mean.human = (sum(callBS == "True")+1)/(n()+2),
            se.human = sqrt((mean.human*(1-mean.human)/n())),
            n=n()) %>%
  mutate(role = "Detector") %>%
  select(role, p, k , mean.human, se.human)
humanComp <- bind_rows(humanCompLie, humanCompDetect)

modelVsHuman <- left_join(modelComp, humanComp, by=c("role", "p", "k"))

comp <- modelVsHuman %>%
  mutate(p.txt = paste("p =", p)) %>%
  ggplot(aes(x=mean.model, y=mean.human, colour=p.txt, shape=role)) +
  geom_point() +
  geom_errorbar(aes(ymin=mean.human-se.human, ymax=mean.human+se.human), width=.02) +
  #geom_errorbarh(aes(xmin=mean.model-se.model, xmax=mean.model+se.model), height=.02) +
  geom_abline(intercept = 0, slope = 1, linetype=2) +
  scale_x_continuous("Model", limits=c(0,1)) +
  scale_y_continuous("Human", limits=c(0,1)) +
  scale_colour_manual(values = condColors.diff) +
  guides(colour=guide_legend(title="")) +
  facet_grid(role~.) +
  theme_minimal() +
  theme(axis.title = element_text(size=18),
        axis.text = element_text(size=18))
comp

df.unif %>%
  filter(role=="Liar", mor==lieMoral, expt=="expt5") %>%
  mutate(p = paste("p =", p)) %>%
  ggplot(aes(x=ks, y=val, colour=p)) +
  geom_abline(size=ablineSize, intercept = 0, slope = 1, colour="darkgray", linetype=2) +
  geom_line(size=lineSize, alpha=alphaLine) +
  ggtitle("Recursive") +
  scale_x_continuous("", labels=rep("",6), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",6), limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))

ggsave("img/modelRecurse_lie_expt5.png", height=2.5, width=3)

df.unif %>%
  filter(role=="Detector", mor==detectMoral, expt=="expt5") %>%
  mutate(p = paste("p =", p)) %>%
  ggplot(aes(x=ks, y=val, colour=p)) +
  geom_line(size=lineSize, alpha=alphaLine) +
  scale_x_continuous("", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))

ggsave("img/modelRecurse_detect_expt5.png", height=2.5, width=3)
df.unif %>%
  filter(role=="Detector", mor==detectMoral, expt=="expt4") %>%
  mutate(p = paste("p =", p)) %>%
  ggplot(aes(x=ks, y=val, colour=p)) +
  geom_line(size=lineSize, alpha=alphaLine) +
  scale_x_continuous("", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))

ggsave("img/modelRecurse_detect_expt4.png", height=2.5, width=3)
df.unif %>%
  filter(role=="Liar", mor==0, expt=="expt4", p==0.5) %>%
  mutate(p = paste("p =", p)) %>%
  ggplot(aes(x=ks, y=val)) +
  geom_line(colour="purple3", size=5) +
  geom_abline(size=4, intercept = 0, slope = 1, colour="darkgray", linetype=2, lwd=1) +
  scale_x_continuous("", breaks=c(0,2,4,6,8,10)) +
  scale_y_continuous("", limits=c(0,10), breaks=c(0,2,4,6,8,10)) +
  theme_minimal() +
  theme(axis.text = element_text(size=24))
## Warning: Duplicated aesthetics after name standardisation: size

ggsave("img/modelRecurse_lie_expt1.png", height=5, width=7)

df.unif %>%
  filter(role=="Detector", mor==0, expt=="expt4", p==0.5) %>%
  mutate(p = paste("p =", p)) %>%
  ggplot(aes(x=ks, y=val)) +
  geom_line(colour="purple3", size=5) +
  scale_x_continuous("", breaks=c(0,2,4,6,8,10)) +
  scale_y_continuous("", limits=c(0,1)) +
  theme_minimal() +
  theme(axis.text = element_text(size=24))

ggsave("img/modelRecurse_detect_expt1.png", height=5, width=7)

MLE Analysis

set.seed(100)
logitToProb <- function(logit){
  exp(logit) / (1+exp(logit))
}

probToLogit <- function(prob){
  log(prob / (1 - prob))
}

quadratic <- function(x, a, b, c){
  a*(x-b)^2+c
}

logisticModel <- function(x, a, b, c){
  logitToProb(pmin(10, pmax(-10, quadratic(x, a, b, c))))
}

lapsePlusLogistic <- function(x, a, b, c, alph){
  alph2 = logitToProb(alph)
  a = exp(a)
  alph2/2 + (1-alph2)*logisticModel(x, a, b, c)
}

loglik <- function(x, y, a, b, c, alph){
  sum(
    dbinom(y, 1, lapsePlusLogistic(x, a, b, c, alph), log=T)
  )
}

brutefit <- function(tmp){
  nLL <- function(a, b, c, alph){
    -loglik(tmp$reportedDrawn, tmp$callBS, a, b, c, alph)+
      (a-1)^2+
      (b-5)^2+
      c^2+
      (alph+2)^2
  }

  iter = 0
  fits = NULL
  fit = NULL
  while(is.null(fits)){
    try(fit <- summary(mle(nLL,
                           start=list(a=rnorm(1, 0, 3),
                                      b=rnorm(1, 0, 3),
                                      c=rnorm(1, 0, 3),
                                      alph=rnorm(1, -2, 3)), method = 'BFGS'), TRUE))
    iter = iter+1

    if(! is.null(fit)){
      fits <- c(tmp$probabilityRed[1], -0.5*fit@m2logL, length(tmp$reportedDrawn), fit@coef[,"Estimate"], fit@coef[,"Std. Error"])
    } else {
      if(iter>100){
        fits <- c(tmp$probabilityRed[1], -9999, 0, 0, 0, 0, 0, 0, 0, 0, 0)
      }
    }
  }
  names(fits) <- c("probabilityRed", "logL", "n", "est.a", "est.b", "est.c", "est.alph", "se.a", "se.b", "se.c", "se.alph")
  return(fits)
}

ps = c(-0.5, 0.5, -0.5, 0.5, -0.5, 0.5, -5, -0.5)
names(ps) <- c("min_a", "max_a", "min_b", "max_b", "min_c", "max_c", "min_alph", "max_alph")

humanDetect_expt4 <- bs.final %>%
  filter(roleCurrent == "bullshitDetector", expt=="expt4") %>%
  mutate(callBS = ifelse(callBS == "False", FALSE, TRUE),
         avgRed = probabilityRed * 10,
         probabilityRed.f = relevel(as.factor(probabilityRed), ref="0.5"))
quadr.fits_expt4 <- data.frame(do.call(rbind, by(humanDetect_expt4, humanDetect_expt4$probabilityRed.f, brutefit)))
print(paste("Failed quadratic fits:", sum(quadr.fits_expt4$logL==-9999)))
## [1] "Failed quadratic fits: 0"
quadr.est_expt4 <- quadr.fits_expt4 %>%
  select(1:7) %>%
  gather("variable", "estimate", 4:7) %>%
  mutate(variable = gsub("est.", "", variable))
quadr.se_expt4 <- quadr.fits_expt4 %>%
  select(-c(4:7)) %>%
  gather("variable", "std.err", 4:7) %>%
  mutate(variable = gsub("se.", "", variable))
quadr.summ_expt4 <- left_join(quadr.est_expt4, quadr.se_expt4, by=c("probabilityRed","logL","n","variable"))
quadr.summ_expt4 %>%
  mutate(z.value = estimate/std.err,
         p.value = ifelse(pnorm(-abs(z.value)) > .000001, pnorm(-abs(z.value)), "<.000001")) %>%
  arrange(probabilityRed, variable)
##    probabilityRed       logL    n variable   estimate   std.err    z.value
## 1             0.2  -759.3295 1450        a -1.5787144 0.2162539  -7.300283
## 2             0.2  -759.3295 1450     alph -1.0657852 0.2254409  -4.727560
## 3             0.2  -759.3295 1450        b  2.4973916 0.1780246  14.028349
## 4             0.2  -759.3295 1450        c -2.4566782 0.3115967  -7.884161
## 5             0.5  -978.7892 1650        a -1.3757927 0.1567622  -8.776305
## 6             0.5  -978.7892 1650     alph -0.6656776 0.1880084  -3.540681
## 7             0.5  -978.7892 1650        b  3.6201312 0.1034947  34.978899
## 8             0.5  -978.7892 1650        c -2.8164056 0.3744991  -7.520460
## 9             0.8 -1077.7395 1650        a -2.4318840 0.1773933 -13.708993
## 10            0.8 -1077.7395 1650     alph -1.5959667 0.7590546  -2.102572
## 11            0.8 -1077.7395 1650        b  4.6323730 0.1488654  31.117853
## 12            0.8 -1077.7395 1650        c -1.5768283 0.2778555  -5.674994
##                 p.value
## 1              <.000001
## 2  1.13616886288692e-06
## 3              <.000001
## 4              <.000001
## 5              <.000001
## 6  0.000199548057724231
## 7              <.000001
## 8              <.000001
## 9              <.000001
## 10     0.01775161530653
## 11             <.000001
## 12             <.000001
# y = a(x-b)^2 + c
b_0.2_expt4 <- quadr.summ_expt4 %>%
  filter(variable == "b", probabilityRed == 0.2)
b_0.5_expt4 <- quadr.summ_expt4 %>%
  filter(variable == "b", probabilityRed == 0.5)
b_0.8_expt4 <- quadr.summ_expt4 %>%
  filter(variable == "b", probabilityRed == 0.8)

# one sided 2 sample t test with pooled variance
wald.z.test <- function(m1, sd1, m2, sd2){
  z <- (m1 - m2) / sqrt(sd1^2 + sd2^2)
  p <- pnorm(abs(z), lower.tail=F)
  return(data.frame(m1, sd1, m2, sd2, z, p))
}

wald.z.test(b_0.5_expt4$estimate, b_0.5_expt4$std.err, b_0.2_expt4$estimate, b_0.2_expt4$std.err)
##         m1       sd1       m2       sd2        z            p
## 1 3.620131 0.1034947 2.497392 0.1780246 5.452254 2.486771e-08
wald.z.test(b_0.5_expt4$estimate, b_0.5_expt4$std.err, b_0.8_expt4$estimate, b_0.8_expt4$std.err)
##         m1       sd1       m2       sd2         z            p
## 1 3.620131 0.1034947 4.632373 0.1488654 -5.583039 1.181758e-08
quadr.plot_expt4 <- data.frame(x=rep(0:10,3), 
           p=as.factor(c(rep(0.2,11), rep(0.5,11), rep(0.8,11))),
           a=c(rep(quadr.fits_expt4$est.a[quadr.fits_expt4$probabilityRed==0.2],11),
               rep(quadr.fits_expt4$est.a[quadr.fits_expt4$probabilityRed==0.5],11),
               rep(quadr.fits_expt4$est.a[quadr.fits_expt4$probabilityRed==0.8],11)),
           b=c(rep(quadr.fits_expt4$est.b[quadr.fits_expt4$probabilityRed==0.2],11),
               rep(quadr.fits_expt4$est.b[quadr.fits_expt4$probabilityRed==0.5],11),
               rep(quadr.fits_expt4$est.b[quadr.fits_expt4$probabilityRed==0.8],11)),
           c=c(rep(quadr.fits_expt4$est.c[quadr.fits_expt4$probabilityRed==0.2],11),
               rep(quadr.fits_expt4$est.c[quadr.fits_expt4$probabilityRed==0.5],11),
               rep(quadr.fits_expt4$est.c[quadr.fits_expt4$probabilityRed==0.8],11)),
           alph=c(rep(quadr.fits_expt4$est.alph[quadr.fits_expt4$probabilityRed==0.2],11),
               rep(quadr.fits_expt4$est.alph[quadr.fits_expt4$probabilityRed==0.5],11),
               rep(quadr.fits_expt4$est.alph[quadr.fits_expt4$probabilityRed==0.8],11))) %>%
  mutate(prob=lapsePlusLogistic(x, a, b, c, alph))
ggplot(quadr.plot_expt4, aes(x=x, y=prob, colour=p)) +
  geom_line() +
  scale_y_continuous(limits=c(0,1))

ggsave("img/quadraticFit_expt4.png")
## Saving 7 x 5 in image
exptDetect_expt4 <- bs.final %>%
  filter(expt=="expt4", roleCurrent == "bullshitDetector") %>%
  mutate(probabilityRed.txt = paste("p =", probabilityRed)) %>%
  group_by(probabilityRed.txt, reportedDrawn) %>%
  summarise(prop = (sum(callBS == "True")+1)/(n()+2),
            se = sqrt((prop*(1-prop)/n())),
            n = n()) %>%
  filter(n > 3) %>%
  mutate(fit="experiment")
quadr.plot2_expt4 <- data.frame(probabilityRed.txt = paste("p =", quadr.plot_expt4$p),
                          reportedDrawn = quadr.plot_expt4$x,
                          prop = quadr.plot_expt4$prob,
                          se = NA,
                          n = NA,
                          fit = "model")
quadr.plot_full_expt4 <- bind_rows(exptDetect_expt4, quadr.plot2_expt4)
ggplot(data=quadr.plot_full_expt4, aes(x=reportedDrawn, y=prop, colour=probabilityRed.txt, linetype=fit)) +
  geom_point(data=filter(quadr.plot_full_expt4, fit=="experiment")) +
  geom_line() +
  geom_errorbar(data=filter(quadr.plot_full_expt4, fit=="experiment"), aes(x=reportedDrawn, colour=probabilityRed.txt, min=prop-se, max=prop+se), width=.3) +
  scale_x_continuous("Reported Marbles Drawn", limits=c(-0.2,10.2)) +
  scale_y_continuous("Proportion BS Called", limits=c(0,1)) +
  guides(colour=guide_legend(title="")) +
  theme_minimal()

ggsave("img/fit_expt_model_expt4.png")
## Saving 7 x 5 in image
humanDetect_expt5 <- bs.final %>%
  filter(roleCurrent == "bullshitDetector", expt=="expt5") %>%
  mutate(callBS = ifelse(callBS == "False", FALSE, TRUE),
         avgRed = probabilityRed * 10,
         probabilityRed.f = relevel(as.factor(probabilityRed), ref="0.5"))
quadr.fits_expt5 <- data.frame(do.call(rbind, by(humanDetect_expt5, humanDetect_expt5$probabilityRed.f, brutefit)))
print(paste("Failed quadratic fits:", sum(quadr.fits_expt5$logL==-9999)))
## [1] "Failed quadratic fits: 0"
quadr.est_expt5 <- quadr.fits_expt5 %>%
  select(1:7) %>%
  gather("variable", "estimate", 4:7) %>%
  mutate(variable = gsub("est.", "", variable))
quadr.se_expt5 <- quadr.fits_expt5 %>%
  select(-c(4:7)) %>%
  gather("variable", "std.err", 4:7) %>%
  mutate(variable = gsub("se.", "", variable))
quadr.summ_expt5 <- left_join(quadr.est_expt5, quadr.se_expt5, by=c("probabilityRed","logL","n","variable"))
quadr.summ_expt5 %>%
  mutate(z.value = estimate/std.err,
         p.value = ifelse(pnorm(-abs(z.value)) > .000001, pnorm(-abs(z.value)), "<.000001")) %>%
  arrange(probabilityRed, variable)
##    probabilityRed       logL    n variable   estimate   std.err    z.value
## 1             0.2 -1077.2295 1650        a -2.4779787 0.1414685 -17.516113
## 2             0.2 -1077.2295 1650     alph -1.8713911 0.6878112  -2.720792
## 3             0.2 -1077.2295 1650        b  5.2304619 0.1459466  35.838186
## 4             0.2 -1077.2295 1650        c -1.5506620 0.2065455  -7.507605
## 5             0.5  -764.6721 1200        a -1.8089359 0.2464863  -7.338889
## 6             0.5  -764.6721 1200     alph -0.6161343 0.3583381  -1.719422
## 7             0.5  -764.6721 1200        b  6.2511026 0.1675899  37.299992
## 8             0.5  -764.6721 1200        c -1.8466357 0.3751881  -4.921894
## 9             0.8  -935.9940 1600        a -1.0356535 0.1653662  -6.262787
## 10            0.8  -935.9940 1600     alph -0.5335905 0.1455806  -3.665259
## 11            0.8  -935.9940 1600        b  8.0767320 0.1470767  54.915121
## 12            0.8  -935.9940 1600        c -2.4498410 0.2654709  -9.228285
##                 p.value
## 1              <.000001
## 2   0.00325628799493254
## 3              <.000001
## 4              <.000001
## 5              <.000001
## 6    0.0427687901995322
## 7              <.000001
## 8              <.000001
## 9              <.000001
## 10 0.000123543971499205
## 11             <.000001
## 12             <.000001
# y = a(x-b)^2 + c
b_0.2_expt5 <- quadr.summ_expt5 %>%
  filter(variable == "b", probabilityRed == 0.2)
b_0.5_expt5 <- quadr.summ_expt5 %>%
  filter(variable == "b", probabilityRed == 0.5)
b_0.8_expt5 <- quadr.summ_expt5 %>%
  filter(variable == "b", probabilityRed == 0.8)

wald.z.test(b_0.5_expt5$estimate, b_0.5_expt5$std.err, b_0.2_expt5$estimate, b_0.2_expt5$std.err)
##         m1       sd1       m2       sd2        z            p
## 1 6.251103 0.1675899 5.230462 0.1459466 4.592693 2.187808e-06
wald.z.test(b_0.5_expt5$estimate, b_0.5_expt5$std.err, b_0.8_expt5$estimate, b_0.8_expt5$std.err)
##         m1       sd1       m2       sd2         z            p
## 1 6.251103 0.1675899 8.076732 0.1470767 -8.187591 1.332534e-16
quadr.plot_expt5 <- data.frame(x=rep(0:10,3), 
           p=as.factor(c(rep(0.2,11), rep(0.5,11), rep(0.8,11))),
           a=c(rep(quadr.fits_expt5$est.a[quadr.fits_expt5$probabilityRed==0.2],11),
               rep(quadr.fits_expt5$est.a[quadr.fits_expt5$probabilityRed==0.5],11),
               rep(quadr.fits_expt5$est.a[quadr.fits_expt5$probabilityRed==0.8],11)),
           b=c(rep(quadr.fits_expt5$est.b[quadr.fits_expt5$probabilityRed==0.2],11),
               rep(quadr.fits_expt5$est.b[quadr.fits_expt5$probabilityRed==0.5],11),
               rep(quadr.fits_expt5$est.b[quadr.fits_expt5$probabilityRed==0.8],11)),
           c=c(rep(quadr.fits_expt5$est.c[quadr.fits_expt5$probabilityRed==0.2],11),
               rep(quadr.fits_expt5$est.c[quadr.fits_expt5$probabilityRed==0.5],11),
               rep(quadr.fits_expt5$est.c[quadr.fits_expt5$probabilityRed==0.8],11)),
           alph=c(rep(quadr.fits_expt5$est.alph[quadr.fits_expt5$probabilityRed==0.2],11),
               rep(quadr.fits_expt5$est.alph[quadr.fits_expt5$probabilityRed==0.5],11),
               rep(quadr.fits_expt5$est.alph[quadr.fits_expt5$probabilityRed==0.8],11))) %>%
  mutate(prob=lapsePlusLogistic(x, a, b, c, alph))
ggplot(quadr.plot_expt5, aes(x=x, y=prob, colour=p)) +
  geom_line() +
  scale_y_continuous(limits=c(0,1))

ggsave("img/quadraticFit_expt5.png")
## Saving 7 x 5 in image
exptDetect_expt5 <- bs.final %>%
  filter(expt=="expt5", roleCurrent == "bullshitDetector") %>%
  mutate(probabilityRed.txt = paste("p =", probabilityRed)) %>%
  group_by(probabilityRed.txt, reportedDrawn) %>%
  summarise(prop = (sum(callBS == "True")+1)/(n()+2),
            se = sqrt((prop*(1-prop)/n())),
            n = n()) %>%
  filter(n > 3) %>%
  mutate(fit="experiment")
quadr.plot2_expt5 <- data.frame(probabilityRed.txt = paste("p =", quadr.plot_expt5$p),
                          reportedDrawn = quadr.plot_expt5$x,
                          prop = quadr.plot_expt5$prob,
                          se = NA,
                          n = NA,
                          fit = "model")
quadr.plot_full_expt5 <- bind_rows(exptDetect_expt5, quadr.plot2_expt5)
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
ggplot(data=quadr.plot_full_expt5, aes(x=reportedDrawn, y=prop, colour=probabilityRed.txt, linetype=fit)) +
  geom_point(data=filter(quadr.plot_full_expt5, fit=="experiment")) +
  geom_line() +
  geom_errorbar(data=filter(quadr.plot_full_expt5, fit=="experiment"), aes(x=reportedDrawn, colour=probabilityRed.txt, min=prop-se, max=prop+se), width=.3) +
  scale_x_continuous("Reported Marbles Drawn", limits=c(-0.2,10.2)) +
  scale_y_continuous("Proportion BS Called", limits=c(0,1)) +
  guides(colour=guide_legend(title="")) +
  theme_minimal()

ggsave("img/fit_expt_model_expt5.png")
## Saving 7 x 5 in image
meancomb <- function(mus, sds){
  sum(mus/sds^2)/sum(1/sds^2)
}

## compare utility conditions
chisqu_util <- function(bs, sds){
  b_0.2 = meancomb(c(bs[1],bs[4]), c(sds[1],sds[4]))
  b_0.5 = meancomb(c(bs[2],bs[5]), c(sds[2],sds[5]))
  b_0.8 = meancomb(c(bs[3],bs[6]), c(sds[3],sds[6]))
  sum(((c(bs[1],bs[4]) - b_0.2)/c(sds[1],sds[4]))^2) + sum(((c(bs[2],bs[5]) - b_0.5)/c(sds[2],sds[5]))^2) + sum(((c(bs[3],bs[6]) - b_0.8)/c(sds[3],sds[6]))^2)
}

detectMLE.pred <- data.frame(p=as.factor(rep(c(0.2,0.5,0.8),2)),
                          expt=as.factor(c(rep("expt4",3), rep("expt5",3))),
                          b=c(filter(quadr.summ_expt4, probabilityRed==0.2 & variable=="b")$estimate,
                               filter(quadr.summ_expt4, probabilityRed==0.5 & variable=="b")$estimate,
                               filter(quadr.summ_expt4, probabilityRed==0.8 & variable=="b")$estimate,
                               filter(quadr.summ_expt5, probabilityRed==0.2 & variable=="b")$estimate,
                               filter(quadr.summ_expt5, probabilityRed==0.5 & variable=="b")$estimate,
                               filter(quadr.summ_expt5, probabilityRed==0.8 & variable=="b")$estimate),
                          b.se=c(filter(quadr.summ_expt4, probabilityRed==0.2 & variable=="b")$std.err,
                               filter(quadr.summ_expt4, probabilityRed==0.5 & variable=="b")$std.err,
                               filter(quadr.summ_expt4, probabilityRed==0.8 & variable=="b")$std.err,
                               filter(quadr.summ_expt5, probabilityRed==0.2 & variable=="b")$std.err,
                               filter(quadr.summ_expt5, probabilityRed==0.5 & variable=="b")$std.err,
                               filter(quadr.summ_expt5, probabilityRed==0.8 & variable=="b")$std.err))


matr_b.se <- matrix(detectMLE.pred$b.se, 3, dimnames=list(c("p=0.2","p=0.5","p=0.8"), c("red","blue")))
matr_b.mean <- matrix(detectMLE.pred$b, 3, dimnames=list(c("p=0.2","p=0.5","p=0.8"), c("red","blue")))


(truechi <- chisqu_util(as.vector(matr_b.mean), as.vector(matr_b.se)))
## [1] 590.2743
(pchi <- pchisq(truechi, df=4, lower.tail=FALSE))
## [1] 1.972671e-126
# do thissss
## simulate chi-square
n = 10000
simChiSq <- data.frame(n=numeric(),
                       i=numeric(),
                       j=numeric(),
                       b=numeric(),
                       chisq=numeric())
for(i in 1:n){
  means = rnorm(6,0,as.vector(matr_b.se))
  chi = chisqu_util(means, as.vector(matr_b.se))
  newdf = data.frame(n=i,
                     i=rep(1:3, 2),
                     j=rep(1:2, each=3),
                     mu = means,
                     chisq = chi)
  simChiSq = bind_rows(simChiSq, newdf)
}

orderedSim <- simChiSq %>%
  select(n, chisq) %>%
  group_by(n) %>%
  unique() %>%
  arrange(chisq)

ggplot(orderedSim, aes(x=chisq)) +
  geom_density() +
  geom_vline(xintercept=orderedSim$chisq[9500], colour="red") +
  ggtitle("Simulated chi-square (df=4) distribution")

humanLie <- bs.final %>%
  filter(roleCurrent == "bullshitter") %>%
  mutate(condition = paste0(probabilityRed,"_",expt))

logitToProb <- function(logit){
  exp(logit) / (1+exp(logit))
}

probToLogit <- function(prob){
  log(prob / (1 - prob))
}

# w = b + a
# u = a / (b + a)


lieLogisticModel <- function(k, beta, mu){
  logodds = beta * (k-mu)
  logitToProb(pmin(10, pmax(-10, logodds)))
}

p.lie.k <- function(k, beta, mu){
  lieLogisticModel(k, beta, mu)
}


p.kstar.k <- function(k, kstar, alph){
  alph = logitToProb(alph)
  dbinom(kstar, 10, alph) / (1-dbinom(k, 10, alph))  
}

lieLogisticBinom <- function(k, kstar, beta, mu, alph){
  p.lie = p.lie.k(k, beta, mu)
  ifelse(k == kstar, 1 - p.lie, p.lie * p.kstar.k(k, kstar, alph))
}


#lieLogisticBinom(5,6,2,1,5,2,5,0)

nLL <- function(beta, mu0.2_4, alph0.2_4,
                      mu0.5_4, alph0.5_4,
                      mu0.8_4, alph0.8_4,
                      mu0.2_5, alph0.2_5,
                      mu0.5_5, alph0.5_5,
                      mu0.8_5, alph0.8_5){
  k = humanLie$drawnRed
  kstar = humanLie$reportedDrawn
  expt = humanLie$expt
  prob = humanLie$probabilityRed
  mu = case_when(
    prob == 0.2 & expt == "expt4" ~ mu0.2_4,
    prob == 0.5 & expt == "expt4" ~ mu0.5_4,
    prob == 0.8 & expt == "expt4" ~ mu0.8_4,
    prob == 0.2 & expt == "expt5" ~ mu0.2_5,
    prob == 0.5 & expt == "expt5" ~ mu0.5_5,
    prob == 0.8 & expt == "expt5" ~ mu0.8_5
   )
  alph = case_when(
    prob == 0.2 & expt == "expt4" ~ alph0.2_4,
    prob == 0.5 & expt == "expt4" ~ alph0.5_4,
    prob == 0.8 & expt == "expt4" ~ alph0.8_4,
    prob == 0.2 & expt == "expt5" ~ alph0.2_5,
    prob == 0.5 & expt == "expt5" ~ alph0.5_5,
    prob == 0.8 & expt == "expt5" ~ alph0.8_5
   )
  
  betas = ifelse(expt=="expt4", 1*beta, -1*beta)
  
  pred = lieLogisticBinom(k, kstar, betas, mu, alph)
  # likelihood of observed kstar for that k, given parameters
  neg.log.lik = -1*sum(log(pred))
  mus = c(mu0.2_4, mu0.5_4, mu0.8_4, mu0.2_5, mu0.5_5, mu0.8_5)
  alphas = c(alph0.2_4, alph0.5_4, alph0.8_4, alph0.2_5, alph0.5_5, alph0.8_5)
  #neg.log.prior = sum(.0001*(mus-5)^2)-abs(beta)+ sum(alphas^2)+u*10
  neg.log.prior = 0
  neg.log.lik+neg.log.prior
}

# as bernoulli instead of if statement


fit <- summary(mle(nLL,
           start=list(beta=rnorm(1, 0, 0.5),
                      mu0.2_4=rnorm(1, 5, 3),
                      alph0.2_4=rnorm(1, 0, 0.5),
                      mu0.5_4=rnorm(1, 5, 3),
                      alph0.5_4=rnorm(1, 0, 0.5),
                      mu0.8_4=rnorm(1, 5, 3),
                      alph0.8_4=rnorm(1, 0, 0.5),
                      mu0.2_5=rnorm(1, 5, 3),
                      alph0.2_5=rnorm(1, 0, 0.5),
                      mu0.5_5=rnorm(1, 5, 3),
                      alph0.5_5=rnorm(1, 0, 0.5),
                      mu0.8_5=rnorm(1, 5, 3),
                      alph0.8_5=rnorm(1, 0, 0.5)),
           method = "BFGS"))
fit
## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = nLL, start = list(beta = rnorm(1, 0, 0.5), mu0.2_4 = rnorm(1, 
##     5, 3), alph0.2_4 = rnorm(1, 0, 0.5), mu0.5_4 = rnorm(1, 5, 
##     3), alph0.5_4 = rnorm(1, 0, 0.5), mu0.8_4 = rnorm(1, 5, 3), 
##     alph0.8_4 = rnorm(1, 0, 0.5), mu0.2_5 = rnorm(1, 5, 3), alph0.2_5 = rnorm(1, 
##         0, 0.5), mu0.5_5 = rnorm(1, 5, 3), alph0.5_5 = rnorm(1, 
##         0, 0.5), mu0.8_5 = rnorm(1, 5, 3), alph0.8_5 = rnorm(1, 
##         0, 0.5)), method = "BFGS")
## 
## Coefficients:
##              Estimate Std. Error
## beta      -0.16881301 0.01074292
## mu0.2_4    1.06928791 0.33013399
## alph0.2_4 -0.60546307 0.02646480
## mu0.5_4   -0.53931421 0.46679991
## alph0.5_4  0.11374033 0.02839969
## mu0.8_4    0.25185704 0.55486313
## alph0.8_4  0.75871137 0.03384142
## mu0.2_5    6.43992378 0.38842659
## alph0.2_5 -0.77714666 0.02794583
## mu0.5_5    9.13760619 0.44637813
## alph0.5_5  0.05458649 0.03033608
## mu0.8_5    9.90510701 0.34089381
## alph0.8_5  0.41985170 0.02600065
## 
## -2 log L: 24509.8
binom.p <- function(kstar, alph){
  alph = logitToProb(alph)
  dbinom(kstar, 10, alph)
}

lieMLE.pred <- data.frame(k=rep(rep(0:10,each=11),3*2), 
                          kstar=rep(0:10, 11*3*2),
                          p=as.factor(rep(c(rep(0.2,11*11), rep(0.5,11*11), rep(0.8,11*11)),2)),
                          expt=as.factor(c(rep("expt4",11*11*3), rep("expt5",11*11*3))),
                          betas=c(rep(fit@coef["beta","Estimate"],11*11*3*2)),
                          mu=c(rep(fit@coef["mu0.2_4","Estimate"],11*11),
                               rep(fit@coef["mu0.5_4","Estimate"],11*11),
                               rep(fit@coef["mu0.8_4","Estimate"],11*11),
                               rep(fit@coef["mu0.2_5","Estimate"],11*11),
                               rep(fit@coef["mu0.5_5","Estimate"],11*11),
                               rep(fit@coef["mu0.8_5","Estimate"],11*11)),
                          mu.se=c(rep(fit@coef["mu0.2_4","Std. Error"],11*11),
                                  rep(fit@coef["mu0.5_4","Std. Error"],11*11),
                                  rep(fit@coef["mu0.8_4","Std. Error"],11*11),
                                  rep(fit@coef["mu0.2_5","Std. Error"],11*11),
                                  rep(fit@coef["mu0.5_5","Std. Error"],11*11),
                                  rep(fit@coef["mu0.8_5","Std. Error"],11*11)),
                          alph=c(rep(fit@coef["alph0.2_4","Estimate"],11*11),
                                 rep(fit@coef["alph0.5_4","Estimate"],11*11),
                                 rep(fit@coef["alph0.8_4","Estimate"],11*11),
                                 rep(fit@coef["alph0.2_5","Estimate"],11*11),
                                 rep(fit@coef["alph0.5_5","Estimate"],11*11),
                                 rep(fit@coef["alph0.8_5","Estimate"],11*11)),
                          alph.se=c(rep(fit@coef["alph0.2_4","Std. Error"],11*11),
                                    rep(fit@coef["alph0.5_4","Std. Error"],11*11),
                                    rep(fit@coef["alph0.8_4","Std. Error"],11*11),
                                    rep(fit@coef["alph0.2_5","Std. Error"],11*11),
                                    rep(fit@coef["alph0.5_5","Std. Error"],11*11),
                                    rep(fit@coef["alph0.8_5","Std. Error"],11*11))) %>%
  mutate(betas = ifelse(expt=="expt5", -betas, betas),
         p.lie.k = p.lie.k(k, betas, mu),
         p.kstar.k=lieLogisticBinom(k, kstar, betas, mu, alph),
         binom.kstar = binom.p(kstar, alph),
         logp.kstar.k = log(p.kstar.k))

lieMLE.pred %>%
  group_by(p, expt, k) %>%
  summarise(mean.kstar = sum(kstar*p.kstar.k)) %>%
  ggplot(aes(x=k, y=mean.kstar, colour=p)) +
  geom_line() +
  facet_wrap(~expt)

exptLie <- humanLie %>%
  mutate(probabilityRed.txt = paste("p =", probabilityRed)) %>%
  group_by(probabilityRed.txt, expt, drawnRed, reportedDrawn) %>%
  summarise(n = n(),
            truth = sum(drawnRed == reportedDrawn)) %>%
  ungroup() %>%
  group_by(probabilityRed.txt, expt, drawnRed) %>%
  mutate(p.kstar.k = n/sum(n),
         fit="experiment")
lieMLE.2 <- data.frame(probabilityRed.txt = paste("p =", lieMLE.pred$p),
                               expt = lieMLE.pred$expt,
                               drawnRed = lieMLE.pred$k,
                               reportedDrawn = lieMLE.pred$kstar,
                               p.kstar.k = lieMLE.pred$p.kstar.k,
                               n = NA,
                               fit = "model")
lie.full <- bind_rows(exptLie, lieMLE.2)
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
lie.full %>%
  group_by(fit, probabilityRed.txt, expt, drawnRed) %>%
  summarise(mean.reported = sum(reportedDrawn*p.kstar.k)) %>%
  ggplot(aes(x=drawnRed, y=mean.reported, colour=probabilityRed.txt, linetype=fit)) +
  geom_line() +
  geom_abline(intercept = 0, slope = 1, colour="darkgray", linetype=4) +
  facet_wrap(~expt)

ggsave("img/fit_expt_modelLie.png")
## Saving 7 x 5 in image
lieMLE.3 <- data.frame(probabilityRed.txt = paste("p =", lieMLE.pred$p),
                       expt = lieMLE.pred$expt,
                       drawnRed = lieMLE.pred$k,
                       p.lie.k = lieMLE.pred$p.lie.k,
                       se = NA,
                       n = NA,
                       fit = "model") %>%
  unique()
exptLie.3 <- humanLie %>%
  mutate(probabilityRed.txt = paste("p =", probabilityRed)) %>%
  group_by(probabilityRed.txt, expt, drawnRed) %>%
  summarise(p.lie.k = sum(drawnRed != reportedDrawn)/n(),
            se = sqrt(p.lie.k*(1-p.lie.k)/n()),
            n = n()) %>%
  mutate(fit="experiment")

lie.3.full <- bind_rows(lieMLE.3, exptLie.3)
## Warning in bind_rows_(x, .id): binding factor and character vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): binding factor and character vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
ggplot(data=lie.3.full, aes(x=drawnRed, y=p.lie.k, colour=probabilityRed.txt, linetype=fit)) +
  geom_point(data=filter(lie.3.full, fit=="experiment")) +
  geom_line() +
  geom_errorbar(data=filter(lie.3.full, fit=="experiment"), aes(x=drawnRed, colour=probabilityRed.txt, min=p.lie.k-se, max=p.lie.k+se), width=.3) +
  scale_x_continuous("Actual Marbles Drawn", limits=c(-0.2,10.2)) +
  scale_y_continuous("Proportion Lie", limits=c(0,1)) +
  guides(colour=guide_legend(title="")) +
  facet_wrap(~expt) +
  theme_minimal()

ggsave("img/fit_expt_modelWhenLie.png")
## Saving 7 x 5 in image
x = -20:35
data.frame(p=as.factor(rep(c(rep(0.2,length(x)), rep(0.5,length(x)), rep(0.8,length(x))),2)),
           expt=as.factor(c(rep("expt4",length(x)*3), rep("expt5",length(x)*3))),
           x=rep(x,6),
           beta=fit@coef["beta","Estimate"],
           mu=c(rep(fit@coef["mu0.2_4","Estimate"],length(x)),
                rep(fit@coef["mu0.5_4","Estimate"],length(x)),
                rep(fit@coef["mu0.8_4","Estimate"],length(x)),
                rep(fit@coef["mu0.2_5","Estimate"],length(x)),
                rep(fit@coef["mu0.5_5","Estimate"],length(x)),
                rep(fit@coef["mu0.8_5","Estimate"],length(x)))) %>%
  mutate(beta = ifelse(expt=="expt5", -beta, beta),
         p.lie.k = p.lie.k(x, beta, mu)) %>%
  ggplot(aes(x=x, y=p.lie.k, colour=p)) +
  geom_line() +
  geom_vline(xintercept=0, linetype=2) +
  geom_vline(xintercept=10, linetype=2) +
  scale_y_continuous(limits=c(0,1)) +
  facet_wrap(~expt) +
  theme_minimal()

ggsave("img/fit_expt_modelWhenLieWide.png")
## Saving 7 x 5 in image
lieColors = c("truth" = "springgreen3",
              "lie" = "red3")

propLieGiven <- bs.final %>%
  filter(expt %in% c("expt4","expt5"), roleCurrent == "bullshitter") %>%
  group_by(expt, probabilityRed, drawnRed, reportedDrawn) %>%
  summarise(n=n()) %>%
  ungroup() %>%
  group_by(expt, probabilityRed, drawnRed) %>%
  mutate(total = sum(n),
         prop = n/sum(n),
         meanLie = sum(reportedDrawn * prop),
         logprop = log(n),
         logitprop = log(prop/(1-prop)),
         lie = ifelse(drawnRed != reportedDrawn, "lie", "truth"))

propLieGiven %>%
  filter(expt == "expt4") %>%
  ggplot(aes(x=reportedDrawn, y=prop, fill=lie)) +
  geom_bar(stat="identity") +
  ggtitle("Expt 4 k* | k") +
  scale_fill_manual(values=lieColors) +
  coord_flip() +
  facet_grid(probabilityRed ~ drawnRed)

propLieGiven %>%
  filter(expt == "expt5") %>%
  ggplot(aes(x=reportedDrawn, y=prop, fill=lie)) +
  geom_bar(stat="identity") +
  ggtitle("Expt 5 k* | k") +
  scale_fill_manual(values=lieColors) +
  coord_flip() +
  facet_grid(probabilityRed ~ drawnRed)

Binomial distribution of k* for each condition from MLE parameters

lieMLE.pred %>%
  filter(expt=="expt4") %>%
  mutate(lie = ifelse(k != kstar, "lie", "truth")) %>%
  ggplot(aes(x=kstar, y=p.kstar.k, fill=lie)) +
  geom_bar(stat="identity") +
  ggtitle("Expt 4 k* | k") +
  scale_fill_manual(values=lieColors) +
  coord_flip() +
  facet_grid(p ~ k)

ggsave("img/fit_expt_modelLieSpike_expt4.png")
## Saving 7 x 5 in image
lieMLE.pred %>%
  filter(expt=="expt5") %>%
  mutate(lie = ifelse(k != kstar, "lie", "truth")) %>%
  ggplot(aes(x=kstar, y=p.kstar.k, fill=lie)) +
  geom_bar(stat="identity") +
  ggtitle("Expt 5 k* | k") +
  scale_fill_manual(values=lieColors) +
  coord_flip() +
  facet_grid(p ~ k)

ggsave("img/fit_expt_modelLieSpike_expt5.png")
## Saving 7 x 5 in image
ggplot(lieMLE.pred, aes(x=kstar, y=binom.kstar, colour=p)) +
  geom_line() +
  facet_grid(expt~.)

lieMLE.pred %>%
  select(p, expt, alph, alph.se) %>%
  unique() %>%
  mutate(p_expt = paste0(p, "_", expt),
         meanAlph = 10*logitToProb(alph)) %>%
  ggplot(aes(x=p, y=meanAlph, fill=expt)) +
  geom_bar(stat="identity", position="dodge") +
  geom_errorbar(aes(x=p, fill=expt, min=10*logitToProb(alph-alph.se), max=10*logitToProb(alph+alph.se)), width=.3, position=position_dodge(.9)) +
  scale_y_continuous("Mean Lie", limits=c(0,10))
## Warning: Ignoring unknown aesthetics: fill

MLEparams <- lieMLE.pred %>%
  group_by(p, expt) %>%
  select(mu, mu.se, alph, alph.se) %>%
  unique()
## Adding missing grouping variables: `p`, `expt`
MLEparams
## # A tibble: 6 x 6
## # Groups:   p, expt [6]
##   p     expt      mu mu.se    alph alph.se
##   <fct> <fct>  <dbl> <dbl>   <dbl>   <dbl>
## 1 0.2   expt4  1.07  0.330 -0.605   0.0265
## 2 0.5   expt4 -0.539 0.467  0.114   0.0284
## 3 0.8   expt4  0.252 0.555  0.759   0.0338
## 4 0.2   expt5  6.44  0.388 -0.777   0.0279
## 5 0.5   expt5  9.14  0.446  0.0546  0.0303
## 6 0.8   expt5  9.91  0.341  0.420   0.0260
# one sided 2 sample t test with pooled variance
wald.z.test <- function(m1, sd1, m2, sd2){
  z <- (m1 - m2) / sqrt(sd1^2 + sd2^2)
  p <- pnorm(abs(z), lower.tail=F) 
  signif <- p < .05
  return(data.frame(m1, sd1, m2, sd2, z, p, signif))
}

#Comparing mus

# ANOVA

matr_mu.se <- matrix(MLEparams$mu.se, 3, dimnames=list(c("p=0.2","p=0.5","p=0.8"), c("red","blue")))
matr_mu.mean <- matrix(MLEparams$mu, 3, dimnames=list(c("p=0.2","p=0.5","p=0.8"), c("red","blue")))
matr_alph.se <- matrix(MLEparams$alph.se, 3, dimnames=list(c("p=0.2","p=0.5","p=0.8"), c("red","blue")))
matr_alph.mean <- matrix(MLEparams$alph, 3, dimnames=list(c("p=0.2","p=0.5","p=0.8"), c("red","blue")))

# alph

## compare utility conditions - z-test
twosampztest <- function(mu1, sd1, mu2, sd2){
  (mu1 - mu2)/sqrt(sd1^2 + sd2^2)
}



meancomb <- function(mus, sds){
  sum(mus/sds^2)/sum(1/sds^2)
}

meancomb(matr_alph.mean[,'red'], matr_alph.se[,'red'])
## [1] -0.01721598
sumsd <- function(sds){
  sqrt(sum(sds^2)/3^2)
}

# differences between utility conditions
z_mean <- twosampztest(mean(matr_alph.mean[,'red']), sumsd(matr_alph.se[,'red']),
             mean(matr_alph.mean[,'blue']), sumsd(matr_alph.se[,'blue']))
z_mean
## [1] 8.03311
2*pnorm(z_mean, lower.tail=FALSE)
## [1] 9.503234e-16
z_0.2 <- twosampztest(matr_alph.mean['p=0.2','red'], matr_alph.se['p=0.2','red'],
                      matr_alph.mean['p=0.2','blue'], matr_alph.se['p=0.2','blue'])
z_0.2
## [1] 4.460661
2*pnorm(z_0.2, lower.tail=FALSE)
## [1] 8.170736e-06
z_0.5 <- twosampztest(matr_alph.mean['p=0.5','red'], matr_alph.se['p=0.5','red'],
                      matr_alph.mean['p=0.5','blue'], matr_alph.se['p=0.5','blue'])
z_0.5
## [1] 1.423506
2*pnorm(z_0.5, lower.tail=FALSE)
## [1] 0.1545894
z_0.8 <- twosampztest(matr_alph.mean['p=0.8','red'], matr_alph.se['p=0.8','red'],
                      matr_alph.mean['p=0.8','blue'], matr_alph.se['p=0.8','blue'])
z_0.8
## [1] 7.940211
2*pnorm(z_0.8, lower.tail=FALSE)
## [1] 2.018379e-15
## compare probability conditions
chisqu <- function(mus, sds){
  mu_red = meancomb(mus[1:3], sds[1:3])
  mu_blue = meancomb(mus[4:6], sds[4:6])
  sum(((mus[1:3] - mu_red)/sds[1:3])^2) + sum(((mus[4:6] - mu_blue)/sds[4:6])^2)
}

(truechi <- chisqu(as.vector(matr_alph.mean), as.vector(matr_alph.se)))
## [1] 2052.61
(pchi <- pchisq(truechi, df=4, lower.tail=FALSE))
## [1] 0
## simulate chi-square
n = 10000
simChiSq <- data.frame(n=numeric(),
                       i=numeric(),
                       j=numeric(),
                       mu=numeric(),
                       chisq=numeric())
for(i in 1:n){
  means = rnorm(6,0,as.vector(matr_alph.se))
  chi = chisqu(means, as.vector(matr_alph.se))
  newdf = data.frame(n=i,
                     i=rep(1:3, 2),
                     j=rep(1:2, each=3),
                     mu = means,
                     chisq = chi)
  simChiSq = bind_rows(simChiSq, newdf)
}

orderedSim <- simChiSq %>%
  select(n, chisq) %>%
  group_by(n) %>%
  unique() %>%
  arrange(chisq)

ggplot(orderedSim, aes(x=chisq)) +
  geom_density() +
  geom_vline(xintercept=orderedSim$chisq[9500], colour="red") +
  ggtitle("Simulated chi-square (df=4) distribution")

ggplot(orderedSim, aes(x=chisq)) +
  geom_density() +
  geom_vline(xintercept=truechi, linetype=2) +
  scale_x_continuous(limits=c(0,2500))

max(orderedSim$chisq)
## [1] 22.09275
max(orderedSim$chisq) < truechi
## [1] TRUE
# true chi-square distribution
chisample <- data.frame(x=rchisq(100000, 4)) %>%
  arrange(x)
ggplot(chisample, aes(x=x)) +
  geom_density() +
  geom_vline(xintercept=chisample$x[95000], colour="red") +
  ggtitle("Chi-square (df=4) distribution")

# no differences between mus (frequency of lying)?

z_mean_mu <- twosampztest(mean(matr_mu.mean[,'red']), sumsd(matr_mu.se[,'red']),
             mean(matr_mu.mean[,'blue']), sumsd(matr_mu.se[,'blue']))
z_mean_mu
## [1] -23.53952
2*pnorm(abs(z_mean_mu), lower.tail=FALSE)
## [1] 1.60747e-122
(truechi_mu <- chisqu(as.vector(matr_mu.mean), as.vector(matr_mu.se)))
## [1] 55.1496
(pchi_mu <- pchisq(truechi_mu, df=4, lower.tail=FALSE))
## [1] 3.02274e-11
## simulate chi-square
n = 10000
simChiSq_mu <- data.frame(n=numeric(),
                       i=numeric(),
                       j=numeric(),
                       mu=numeric(),
                       chisq=numeric())
for(i in 1:n){
  means = rnorm(6,0,as.vector(matr_mu.se))
  chi = chisqu(means, as.vector(matr_mu.se))
  newdf = data.frame(n=i,
                     i=rep(1:3, 2),
                     j=rep(1:2, each=3),
                     mu = means,
                     chisq = chi)
  simChiSq_mu = bind_rows(simChiSq_mu, newdf)
}

orderedSim_mu <- simChiSq_mu %>%
  select(n, chisq) %>%
  group_by(n) %>%
  unique() %>%
  arrange(chisq)

ggplot(orderedSim_mu, aes(x=chisq)) +
  geom_density() +
  geom_vline(xintercept=orderedSim_mu$chisq[9500], colour="red")

ggplot(orderedSim_mu, aes(x=chisq)) +
  geom_density() +
  geom_vline(xintercept=truechi_mu, linetype=2) +
  scale_x_continuous(limits=c(0,75))

Model predictions about liar’s likelihood to lie

# It is not very intuitive what this means
moral = 0
Expt = 4
recurse_propLie_expt4 <- data.frame(role = "Liar",
                              model="recursive ",
                              p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                              ks=rep(0:10,3),
                              expt = "expt4",
                              val=c(p.lie.k.fromDet(0.2, filter(df.unif, role=="Detector" & mor==0 & expt=="expt4" & p==0.2)$val),
                                    p.lie.k.fromDet(0.5, filter(df.unif, role=="Detector" & mor==0 & expt=="expt4" & p==0.5)$val),
                                    p.lie.k.fromDet(0.8, filter(df.unif, role=="Detector" & mor==0 & expt=="expt4" & p==0.8)$val)))
Expt = 5
recurse_propLie_expt5 <- data.frame(role = "Liar",
                              model="recursive ",
                              p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                              ks=rep(0:10,3),
                              expt = "expt5",
                              val=c(p.lie.k.fromDet(0.2, filter(df.unif, role=="Detector" & mor==0 & expt=="expt5" & p==0.2)$val),
                                    p.lie.k.fromDet(0.5, filter(df.unif, role=="Detector" & mor==0 & expt=="expt5" & p==0.5)$val),
                                    p.lie.k.fromDet(0.8, filter(df.unif, role=="Detector" & mor==0 & expt=="expt5" & p==0.8)$val)))
recurse_propLie <- bind_rows(recurse_propLie_expt4, recurse_propLie_expt5)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
moral = 20
Expt = 4
averse_propLie_expt4 <- data.frame(role = "Liar",
                              model="recursive+aversion",
                              p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                              ks=rep(0:10,3),
                              expt = "expt4",
                              val=c(p.lie.k.fromDet(0.2, filter(df.unif, role=="Detector" & mor==20 & expt=="expt4" & p==0.2)$val),
                                    p.lie.k.fromDet(0.5, filter(df.unif, role=="Detector" & mor==20 & expt=="expt4" & p==0.5)$val),
                                    p.lie.k.fromDet(0.8, filter(df.unif, role=="Detector" & mor==20 & expt=="expt4" & p==0.8)$val)))
Expt = 5
averse_propLie_expt5 <- data.frame(role = "Liar",
                              model="recursive+aversion",
                              p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                              ks=rep(0:10,3),
                              expt = "expt5",
                              val=c(p.lie.k.fromDet(0.2, filter(df.unif, role=="Detector" & mor==20 & expt=="expt5" & p==0.2)$val),
                                    p.lie.k.fromDet(0.5, filter(df.unif, role=="Detector" & mor==20 & expt=="expt5" & p==0.5)$val),
                                    p.lie.k.fromDet(0.8, filter(df.unif, role=="Detector" & mor==20 & expt=="expt5" & p==0.8)$val)))
averse_propLie <- bind_rows(averse_propLie_expt4, averse_propLie_expt5)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
moral = 0
Expt = 4
randOpp_propLie_expt4 <- data.frame(role = "Liar",
                              model = "random opponent",
                              p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                              ks=rep(0:10,3),
                              expt = "expt4",
                              val=c(p.lie.k.fromDet(0.2, rep(bonddepaulo, 11)),
                                    p.lie.k.fromDet(0.5, rep(bonddepaulo, 11)),
                                    p.lie.k.fromDet(0.8, rep(bonddepaulo, 11))))
Expt = 5
randOpp_propLie_expt5 <- data.frame(role = "Liar",
                              model = "random opponent",
                              p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
                              ks=rep(0:10,3),
                              expt = "expt5",
                              val=c(p.lie.k.fromDet(0.2, rep(bonddepaulo, 11)),
                                    p.lie.k.fromDet(0.5, rep(bonddepaulo, 11)),
                                    p.lie.k.fromDet(0.8, rep(bonddepaulo, 11))))
randOpp_propLie <- bind_rows(randOpp_propLie_expt4, randOpp_propLie_expt5)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
propLie.models <- bind_rows(randOpp_propLie, recurse_propLie, averse_propLie) %>%
  mutate(p=paste0("p = ", p))
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
propLies.modelsFig <- propLie.models %>%
  filter(expt == "expt4") %>%
  ggplot(aes(x=ks, y=val, colour=p)) +
  geom_line() +
  scale_x_continuous("Actual Marbles Drawn") +
  scale_y_continuous("Proportion Lie", limits=c(0,1)) +
  scale_colour_manual(values = condColors.diff, guide="none") +
  facet_wrap(~model) +
  theme_bw()

propLies_wFit <- lie.3.full %>%
  filter(expt == "expt4") %>%
  ggplot(aes(x=drawnRed, y=p.lie.k, colour=probabilityRed.txt, linetype=fit)) +
  geom_point(data=filter(lie.3.full, fit=="experiment" & expt == "expt4")) +
  geom_line() +
  geom_errorbar(data=filter(lie.3.full, fit=="experiment" & expt == "expt4"), aes(x=drawnRed, colour=probabilityRed.txt, min=p.lie.k-se, max=p.lie.k+se), width=.3) +
  ggtitle("Human Results") +
  scale_x_continuous("", limits=c(-0.2,10.2)) +
  scale_y_continuous("", limits=c(0,1)) +
  scale_colour_manual(values = condColors.diff, guide="none") +
  guides(colour=guide_legend(title="")) +
  theme_minimal() +
  theme(legend.position = "none")
propLies_wFit

plot_grid(propLies.modelsFig, propLies_wFit, nrow=1, rel_widths=c(6, 4))

ggsave("img/lierate_expt4.png")
## Saving 7 x 5 in image
propLies.modelsFig5 <- propLie.models %>%
  filter(expt == "expt5") %>%
  ggplot(aes(x=ks, y=val, colour=p)) +
  geom_line() +
  scale_x_continuous("Actual Marbles Drawn") +
  scale_y_continuous("Proportion Lie", limits=c(0,1)) +
  scale_colour_manual(values = condColors.diff, guide="none") +
  facet_wrap(~model) +
  theme_bw()

propLies_wFit5 <- lie.3.full %>%
  filter(expt == "expt5") %>%
  ggplot(aes(x=drawnRed, y=p.lie.k, colour=probabilityRed.txt, linetype=fit)) +
  geom_point(data=filter(lie.3.full, fit=="experiment" & expt == "expt5")) +
  geom_line() +
  geom_errorbar(data=filter(lie.3.full, fit=="experiment" & expt == "expt5"), aes(x=drawnRed, colour=probabilityRed.txt, min=p.lie.k-se, max=p.lie.k+se), width=.3) +
  ggtitle("Human Results") +
  scale_x_continuous("", limits=c(-0.2,10.2)) +
  scale_y_continuous("", limits=c(0,1)) +
  scale_colour_manual(values = condColors.diff, guide="none") +
  guides(colour=guide_legend(title="")) +
  theme_minimal() +
  theme(legend.position = "none")
propLies_wFit5

plot_grid(propLies.modelsFig5, propLies_wFit5, nrow=1, rel_widths=c(6, 4))

ggsave("img/lierate_expt5.png")
## Saving 7 x 5 in image

Plot mean lie given (p in Binomial distribution)

moral = 0
matr0.2 <- log(p.L_ksay.k.r(0.2, filter(df.unif, role=="Detector" & mor==0 & expt=="expt4" & p==0.2)$val))
hist3D(z=matr0.2, border="black")

matr0.5 <- log(p.L_ksay.k.r(0.5, filter(df.unif, role=="Detector" & mor==0 & expt=="expt4" & p==0.5)$val))
hist3D(z=matr0.5, border="black")

matr0.8 <- log(p.L_ksay.k.r(0.8, filter(df.unif, role=="Detector" & mor==0 & expt=="expt4" & p==0.8)$val))
hist3D(z=matr0.8, border="black")

#p.L_ksay.k.r(0.5, filter(df.unif, role=="Detector" & mor==0 & expt=="expt4" & p==0.5)$val)

moral = 0
liesFromDetect <- function(model, exptNum, prob){
  if(model=="recurse"){
    detectRate = filter(df.unif, role=="Detector" & mor==0 & expt==paste0("expt",exptNum) & p==prob)$val
  } else if(model=="random opponent"){
    detectRate = rep(bonddepaulo, 11)
  }
  lieMatr <- p.L_ksay.k.r(prob, detectRate)
  data.frame(model=model,
             expt=paste0("expt",exptNum),
             p=prob,
             k=rep(0:10, each=11),
             k.star=rep(0:10, 11),
             prob.k.star = as.vector(lieMatr))
}

Expt = 4
recurseLie_expt4 <- bind_rows(liesFromDetect("recurse", 4, 0.2),
                              liesFromDetect("recurse", 4, 0.5),
                              liesFromDetect("recurse", 4, 0.8))
Expt = 5
recurseLie_expt5 <- bind_rows(liesFromDetect("recurse", 5, 0.2),
                              liesFromDetect("recurse", 5, 0.5),
                              liesFromDetect("recurse", 5, 0.8))
recurseLie <- bind_rows(recurseLie_expt4, recurseLie_expt5) %>%
  mutate(lie = ifelse(k != k.star, "lie", "truth"))

Expt = 4
randomLie_expt4 <- bind_rows(liesFromDetect("random opponent", 4, 0.2),
                             liesFromDetect("random opponent", 4, 0.5),
                             liesFromDetect("random opponent", 4, 0.8))
Expt = 5
randomLie_expt5 <- bind_rows(liesFromDetect("random opponent", 5, 0.2),
                             liesFromDetect("random opponent", 5, 0.5),
                             liesFromDetect("random opponent", 5, 0.8))
randomLie <- bind_rows(randomLie_expt4, randomLie_expt5) %>%
  mutate(lie = ifelse(k != k.star, "lie", "truth"))


recurseLie %>%
  filter(expt == "expt4") %>%
  ggplot(aes(x=k.star, y=prob.k.star, fill=lie)) +
  geom_bar(stat="identity") +
  ggtitle("Recurse Model Expt 4 k* | k") +
  scale_fill_manual(values=lieColors) +
  coord_flip() +
  facet_grid(p ~ k)

recurseLie %>%
  filter(expt == "expt5") %>%
  ggplot(aes(x=k.star, y=prob.k.star, fill=lie)) +
  geom_bar(stat="identity") +
  ggtitle("Recurse Model Expt 5 k* | k") +
  scale_fill_manual(values=lieColors) +
  coord_flip() +
  facet_grid(p ~ k)

randomLie %>%
  filter(expt == "expt4") %>%
  ggplot(aes(x=k.star, y=prob.k.star, fill=lie)) +
  geom_bar(stat="identity") +
  ggtitle("Random Opponent Model Expt 4 k* | k") +
  scale_fill_manual(values=lieColors) +
  coord_flip() +
  facet_grid(p ~ k)

randomLie %>%
  filter(expt == "expt5") %>%
  ggplot(aes(x=k.star, y=prob.k.star, fill=lie)) +
  geom_bar(stat="identity") +
  ggtitle("Random Opponent Model Expt 5 k* | k") +
  scale_fill_manual(values=lieColors) +
  coord_flip() +
  facet_grid(p ~ k)

countRecurseLie <- recurseLie %>%
  mutate(counts = floor(1000*prob.k.star*(0.8*dbinom(k, 10, p))+0.2*(1/11))) %>%
  uncount(weights=counts)

nLL_recurse <- function(beta, mu0.2_4, alph0.2_4,
                      mu0.5_4, alph0.5_4,
                      mu0.8_4, alph0.8_4,
                      mu0.2_5, alph0.2_5,
                      mu0.5_5, alph0.5_5,
                      mu0.8_5, alph0.8_5){
  k = countRecurseLie$k
  kstar = countRecurseLie$k.star
  expt = countRecurseLie$expt
  prob = countRecurseLie$p
  mu = case_when(
    prob == 0.2 & expt == "expt4" ~ mu0.2_4,
    prob == 0.5 & expt == "expt4" ~ mu0.5_4,
    prob == 0.8 & expt == "expt4" ~ mu0.8_4,
    prob == 0.2 & expt == "expt5" ~ mu0.2_5,
    prob == 0.5 & expt == "expt5" ~ mu0.5_5,
    prob == 0.8 & expt == "expt5" ~ mu0.8_5
   )
  alph = case_when(
    prob == 0.2 & expt == "expt4" ~ alph0.2_4,
    prob == 0.5 & expt == "expt4" ~ alph0.5_4,
    prob == 0.8 & expt == "expt4" ~ alph0.8_4,
    prob == 0.2 & expt == "expt5" ~ alph0.2_5,
    prob == 0.5 & expt == "expt5" ~ alph0.5_5,
    prob == 0.8 & expt == "expt5" ~ alph0.8_5
   )
  
  betas = ifelse(expt=="expt4", 1*beta, -1*beta)
  
  pred = lieLogisticBinom(k, kstar, betas, mu, alph)
  # likelihood of observed kstar for that k, given parameters
  neg.log.lik = -1*sum(log(pred))
  mus = c(mu0.2_4, mu0.5_4, mu0.8_4, mu0.2_5, mu0.5_5, mu0.8_5)
  alphas = c(alph0.2_4, alph0.5_4, alph0.8_4, alph0.2_5, alph0.5_5, alph0.8_5)
  #neg.log.prior = sum(.0001*(mus-5)^2)-abs(beta)+ sum(alphas^2)+u*10
  neg.log.prior = 0
  neg.log.lik+neg.log.prior
}



fit_recurse <- summary(mle(nLL_recurse,
           start=list(beta=rnorm(1, 0, 0.5),
                      mu0.2_4=rnorm(1, 5, 3),
                      alph0.2_4=rnorm(1, 0, 0.5),
                      mu0.5_4=rnorm(1, 5, 3),
                      alph0.5_4=rnorm(1, 0, 0.5),
                      mu0.8_4=rnorm(1, 5, 3),
                      alph0.8_4=rnorm(1, 0, 0.5),
                      mu0.2_5=rnorm(1, 5, 3),
                      alph0.2_5=rnorm(1, 0, 0.5),
                      mu0.5_5=rnorm(1, 5, 3),
                      alph0.5_5=rnorm(1, 0, 0.5),
                      mu0.8_5=rnorm(1, 5, 3),
                      alph0.8_5=rnorm(1, 0, 0.5)),
           method = "BFGS"))
fit_recurse
## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = nLL_recurse, start = list(beta = rnorm(1, 0, 
##     0.5), mu0.2_4 = rnorm(1, 5, 3), alph0.2_4 = rnorm(1, 0, 0.5), 
##     mu0.5_4 = rnorm(1, 5, 3), alph0.5_4 = rnorm(1, 0, 0.5), mu0.8_4 = rnorm(1, 
##         5, 3), alph0.8_4 = rnorm(1, 0, 0.5), mu0.2_5 = rnorm(1, 
##         5, 3), alph0.2_5 = rnorm(1, 0, 0.5), mu0.5_5 = rnorm(1, 
##         5, 3), alph0.5_5 = rnorm(1, 0, 0.5), mu0.8_5 = rnorm(1, 
##         5, 3), alph0.8_5 = rnorm(1, 0, 0.5)), method = "BFGS")
## 
## Coefficients:
##             Estimate Std. Error
## beta      -0.6400938 0.02887682
## mu0.2_4    3.8794087 0.15490698
## alph0.2_4  0.1811010 0.02755201
## mu0.5_4    4.8209756 0.12473705
## alph0.5_4  0.2430890 0.03254290
## mu0.8_4    6.0218842 0.15737649
## alph0.8_4  0.9486908 0.04869127
## mu0.2_5    3.7987160 0.15211178
## alph0.2_5 -0.9981020 0.04758842
## mu0.5_5    5.0803793 0.12466849
## alph0.5_5 -0.2937834 0.03242535
## mu0.8_5    6.2380700 0.15086714
## alph0.8_5 -0.1254952 0.02768063
## 
## -2 log L: 17967.7
lieMLE.pred_recurse <- data.frame(k=rep(rep(0:10,each=11),3*2), 
                          kstar=rep(0:10, 11*3*2),
                          p=as.factor(rep(c(rep(0.2,11*11), rep(0.5,11*11), rep(0.8,11*11)),2)),
                          expt=as.factor(c(rep("expt4",11*11*3), rep("expt5",11*11*3))),
                          alph=c(rep(fit_recurse@coef["alph0.2_4","Estimate"],11*11),
                                 rep(fit_recurse@coef["alph0.5_4","Estimate"],11*11),
                                 rep(fit_recurse@coef["alph0.8_4","Estimate"],11*11),
                                 rep(fit_recurse@coef["alph0.2_5","Estimate"],11*11),
                                 rep(fit_recurse@coef["alph0.5_5","Estimate"],11*11),
                                 rep(fit_recurse@coef["alph0.8_5","Estimate"],11*11)),
                          alph.se=c(rep(fit_recurse@coef["alph0.2_4","Std. Error"],11*11),
                                    rep(fit_recurse@coef["alph0.5_4","Std. Error"],11*11),
                                    rep(fit_recurse@coef["alph0.8_4","Std. Error"],11*11),
                                    rep(fit_recurse@coef["alph0.2_5","Std. Error"],11*11),
                                    rep(fit_recurse@coef["alph0.5_5","Std. Error"],11*11),
                                    rep(fit_recurse@coef["alph0.8_5","Std. Error"],11*11))) %>%
  mutate(p.kstar.k = p.kstar.k(k, kstar, alph),
         binom.kstar = binom.p(kstar, alph),
         alph.asProb = logitToProb(alph))

ggplot(lieMLE.pred_recurse, aes(x=kstar, y=binom.kstar, colour=p)) +
  geom_line() +
  facet_grid(expt~.)

countRandomLie <- randomLie %>%
  mutate(counts = floor(1000*prob.k.star*(0.8*dbinom(k, 10, p))+0.2*(1/11))) %>%
  uncount(weights=counts)

nLL_random <- function(beta, mu0.2_4, alph0.2_4,
                      mu0.5_4, alph0.5_4,
                      mu0.8_4, alph0.8_4,
                      mu0.2_5, alph0.2_5,
                      mu0.5_5, alph0.5_5,
                      mu0.8_5, alph0.8_5){
  k = countRandomLie$k
  kstar = countRandomLie$k.star
  expt = countRandomLie$expt
  prob = countRandomLie$p
  mu = case_when(
    prob == 0.2 & expt == "expt4" ~ mu0.2_4,
    prob == 0.5 & expt == "expt4" ~ mu0.5_4,
    prob == 0.8 & expt == "expt4" ~ mu0.8_4,
    prob == 0.2 & expt == "expt5" ~ mu0.2_5,
    prob == 0.5 & expt == "expt5" ~ mu0.5_5,
    prob == 0.8 & expt == "expt5" ~ mu0.8_5
   )
  alph = case_when(
    prob == 0.2 & expt == "expt4" ~ alph0.2_4,
    prob == 0.5 & expt == "expt4" ~ alph0.5_4,
    prob == 0.8 & expt == "expt4" ~ alph0.8_4,
    prob == 0.2 & expt == "expt5" ~ alph0.2_5,
    prob == 0.5 & expt == "expt5" ~ alph0.5_5,
    prob == 0.8 & expt == "expt5" ~ alph0.8_5
   )
  
  betas = ifelse(expt=="expt4", 1*beta, -1*beta)
  
  pred = lieLogisticBinom(k, kstar, betas, mu, alph)
  # likelihood of observed kstar for that k, given parameters
  neg.log.lik = -1*sum(log(pred))
  mus = c(mu0.2_4, mu0.5_4, mu0.8_4, mu0.2_5, mu0.5_5, mu0.8_5)
  alphas = c(alph0.2_4, alph0.5_4, alph0.8_4, alph0.2_5, alph0.5_5, alph0.8_5)
  #neg.log.prior = sum(.0001*(mus-5)^2)-abs(beta)+ sum(alphas^2)+u*10
  neg.log.prior = 0
  neg.log.lik+neg.log.prior
}



fit_random <- summary(mle(nLL_random,
           start=list(beta=rnorm(1, 0, 0.5),
                      mu0.2_4=rnorm(1, 5, 3),
                      alph0.2_4=rnorm(1, 0, 0.5),
                      mu0.5_4=rnorm(1, 5, 3),
                      alph0.5_4=rnorm(1, 0, 0.5),
                      mu0.8_4=rnorm(1, 5, 3),
                      alph0.8_4=rnorm(1, 0, 0.5),
                      mu0.2_5=rnorm(1, 5, 3),
                      alph0.2_5=rnorm(1, 0, 0.5),
                      mu0.5_5=rnorm(1, 5, 3),
                      alph0.5_5=rnorm(1, 0, 0.5),
                      mu0.8_5=rnorm(1, 5, 3),
                      alph0.8_5=rnorm(1, 0, 0.5)),
           method = "BFGS"))
fit_random
## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = nLL_random, start = list(beta = rnorm(1, 0, 0.5), 
##     mu0.2_4 = rnorm(1, 5, 3), alph0.2_4 = rnorm(1, 0, 0.5), mu0.5_4 = rnorm(1, 
##         5, 3), alph0.5_4 = rnorm(1, 0, 0.5), mu0.8_4 = rnorm(1, 
##         5, 3), alph0.8_4 = rnorm(1, 0, 0.5), mu0.2_5 = rnorm(1, 
##         5, 3), alph0.2_5 = rnorm(1, 0, 0.5), mu0.5_5 = rnorm(1, 
##         5, 3), alph0.5_5 = rnorm(1, 0, 0.5), mu0.8_5 = rnorm(1, 
##         5, 3), alph0.8_5 = rnorm(1, 0, 0.5)), method = "BFGS")
## 
## Coefficients:
##             Estimate Std. Error
## beta      -0.5824235 0.03088194
## mu0.2_4    6.4746142 0.30831566
## alph0.2_4  1.1816306 0.02832290
## mu0.5_4    6.7200530 0.16157717
## alph0.5_4  1.1636109 0.03333095
## mu0.8_4    6.8950417 0.14601833
## alph0.8_4  1.1731283 0.04256182
## mu0.2_5    3.1047138 0.14600784
## alph0.2_5 -1.1734258 0.04256525
## mu0.5_5    3.2792532 0.16161302
## alph0.5_5 -1.1638195 0.03333301
## mu0.8_5    3.5176488 0.30913262
## alph0.8_5 -1.1816878 0.02832330
## 
## -2 log L: 19652.7
lieMLE.pred_random <- data.frame(k=rep(rep(0:10,each=11),3*2), 
                          kstar=rep(0:10, 11*3*2),
                          p=as.factor(rep(c(rep(0.2,11*11), rep(0.5,11*11), rep(0.8,11*11)),2)),
                          expt=as.factor(c(rep("expt4",11*11*3), rep("expt5",11*11*3))),
                          alph=c(rep(fit_random@coef["alph0.2_4","Estimate"],11*11),
                                 rep(fit_random@coef["alph0.5_4","Estimate"],11*11),
                                 rep(fit_random@coef["alph0.8_4","Estimate"],11*11),
                                 rep(fit_random@coef["alph0.2_5","Estimate"],11*11),
                                 rep(fit_random@coef["alph0.5_5","Estimate"],11*11),
                                 rep(fit_random@coef["alph0.8_5","Estimate"],11*11)),
                          alph.se=c(rep(fit_random@coef["alph0.2_4","Std. Error"],11*11),
                                    rep(fit_random@coef["alph0.5_4","Std. Error"],11*11),
                                    rep(fit_random@coef["alph0.8_4","Std. Error"],11*11),
                                    rep(fit_random@coef["alph0.2_5","Std. Error"],11*11),
                                    rep(fit_random@coef["alph0.5_5","Std. Error"],11*11),
                                    rep(fit_random@coef["alph0.8_5","Std. Error"],11*11))) %>%
  mutate(p.kstar.k = p.kstar.k(k, kstar, alph),
         binom.kstar = binom.p(kstar, alph),
         alph.asProb = logitToProb(alph))

ggplot(lieMLE.pred_random, aes(x=kstar, y=binom.kstar, colour=p)) +
  geom_line() +
  facet_grid(expt~.)

lieMLE.pred <- lieMLE.pred %>%
  mutate(model="Human")
lieMLE.pred_recurse <- lieMLE.pred_recurse %>%
  mutate(model="Recursive")
lieMLE.pred_random <- lieMLE.pred_random %>%
  mutate(model="Random Opponent")

lieMLE.pred_all <- bind_rows(lieMLE.pred, lieMLE.pred_recurse, lieMLE.pred_random) %>%
  mutate(p = paste("p =", p),
         expt=factor(ifelse(expt=="expt4", "Red", "Blue"), levels=c("Red", "Blue")),
         model=factor(model, levels=c("Random Opponent", "Recursive", "Human")))


lieMLE.pred_human <- lieMLE.pred_all %>%
  filter(model == "Human") %>%
  select(p, expt, alph, alph.se, model) %>%
  unique() %>%
  mutate(meanAlph = 10*logitToProb(alph),
         lowerMeanAlph = 10*logitToProb(alph-2*alph.se),
         upperMeanAlph = 10*logitToProb(alph+2*alph.se))

lieMLE.pred_all %>%
  select(p, expt, alph, alph.se, model) %>%
  unique() %>%
  mutate(meanAlph = 10*logitToProb(alph)) %>%
  ggplot(aes(x=expt, y=meanAlph, fill=p)) +
  geom_bar(stat="identity", position="dodge") +
  geom_errorbar(data=lieMLE.pred_human, aes(x=expt, fill=p, min=lowerMeanAlph, max=upperMeanAlph), width=.3, position=position_dodge(.9)) +
  scale_x_discrete("Utility") +
  scale_y_continuous("Average Lie", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0,0,0.5)) +
  scale_fill_manual(values = condColors.diff) +
  guides(fill=guide_legend(title="", override.aes=list(size=2.5))) +
  facet_grid(model~.) +
  theme_bw()
## Warning: Ignoring unknown aesthetics: fill

ggsave("img/liekparameter.pdf", units="cm", height=13, width=8.7)
liePred_random <- lieMLE.pred_all %>%
  filter(model == "Random Opponent") %>%
  select(p, expt, alph, alph.se) %>%
  unique() %>%
  mutate(meanAlph = 10*logitToProb(alph)) %>%
  ggplot(aes(x=expt, y=meanAlph, fill=p)) +
  geom_bar(stat="identity", position="dodge") +
  ggtitle("Random Opponent") +
  scale_x_discrete("Utility") +
  scale_y_continuous("Average Lie", limits=c(0,10), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_fill_manual(values = condColors.diff) +
  guides(fill=guide_legend(title="", override.aes=list(size=2.5))) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))

liePred_recurse <- lieMLE.pred_all %>%
  filter(model == "Recursive") %>%
  select(p, expt, alph, alph.se) %>%
  unique() %>%
  mutate(meanAlph = 10*logitToProb(alph)) %>%
  ggplot(aes(x=expt, y=meanAlph, fill=p)) +
  geom_bar(stat="identity", position="dodge") +
  ggtitle("Recursive") +
  scale_x_discrete("") +
  scale_y_continuous("", limits=c(0,10), labels=rep("",6), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_fill_manual(values = condColors.diff) +
  guides(fill=guide_legend(title="", override.aes=list(size=2.5))) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))

liePred_human <- lieMLE.pred_all %>%
  filter(model == "Human") %>%
  select(p, expt, alph, alph.se) %>%
  unique() %>%
  mutate(meanAlph = 10*logitToProb(alph)) %>%
  ggplot(aes(x=expt, y=meanAlph, fill=p)) +
  geom_bar(stat="identity", position="dodge") +
  geom_errorbar(data=lieMLE.pred_human, aes(x=expt, fill=p, min=lowerMeanAlph, max=upperMeanAlph), width=.3, position=position_dodge(.9)) +
  ggtitle("Human Results") +
  scale_x_discrete("") +
  scale_y_continuous("", limits=c(0,10), labels=rep("",6), breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_fill_manual(values = condColors.diff) +
  guides(fill=guide_legend(title="", override.aes=list(size=2.5))) +
  theme_minimal() +
  theme(legend.text = element_text(size=7),
        legend.key.size = unit(0.36, 'cm'),
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold"),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
## Warning: Ignoring unknown aesthetics: fill
legend_lie <- get_legend(liePred_human)

liePred <- plot_grid(liePred_random,
                      liePred_recurse,
                      liePred_human + theme(legend.position="none"), 
                      nrow=1,
                      rel_widths=c(10.52, 10, 10))

# liarPredLabel <- ggdraw() +
#   draw_label("Liar", fontface="bold", size=12, y=0.87, x=0.5, hjust=0, angle=270)
liepredrow <- plot_grid(liePred, liarLabel, nrow=1, rel_widths=c(96, 4))

withLegend_lie <- ggdraw() +
  draw_plot(liepredrow) +
  draw_plot(legend_lie, x=0.415, y=0.32)
withLegend_lie

ggsave("img/liePredictions.pdf", withLegend_lie, units="cm", height=5.5, width=17.8)
detect.recurse.nomoral_expt5 <- df.unif %>%
  filter(role=="Detector", mor==0, expt=="expt5") %>%
  mutate(p = paste("p =", p)) %>%
  ggplot(aes(x=ks, y=val, colour=p)) +
  geom_line(size=lineSize, alpha=alphaLine) +
  scale_x_continuous("", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
detect.recurse.nomoral_expt5

ggsave("img/modelRecurse_detect_nomoral_expt5.png", detect.recurse.nomoral_expt5, height=8, width=7.5)


moral = 0
Expt = 5
detect0.2idiot.nomoral_expt5 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.2, lie0.2idiot)
detect0.5idiot.nomoral_expt5 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.5, lie0.5idiot)
detect0.8idiot.nomoral_expt5 <- mapply(p.D_bs.ksay.r, 0:numMarbles, 0.8, lie0.8idiot)

detect.respondToIdiot.nomoral_expt5 <- data.frame(p=c(rep(0.2,11), rep(0.5,11), rep(0.8,11)),
           ksay=rep(0:10,3),
           bs=c(detect0.2idiot_expt5, detect0.5idiot_expt5, detect0.8idiot_expt5)) %>%
  mutate(probabilityRed.txt = paste("p =", p)) %>%
  ggplot(aes(x=ksay, y=bs, colour=probabilityRed.txt)) +
  geom_line(size=lineSize, alpha=alphaLine) +
  # scale_x_continuous("", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  # scale_y_continuous("", labels=rep("",5), limits=c(0,1), expand=c(0,0)) +
  scale_x_continuous("Reported Marbles", breaks=c(0,2,4,6,8,10), expand=c(0,0)) +
  scale_y_continuous("Proportion BS Called", limits=c(0,1), expand=c(0,0)) +
  scale_colour_manual(values=condColors.diff)+
  theme_minimal() +
  theme(legend.position = "none",
        plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize),
        axis.title = element_text(size=axisTitleSize),
        axis.text = element_text(size=axisTextSize),
        axis.text.x = element_text(margin=marginAxisX),
        axis.text.y = element_text(margin=margin(0,-1,0,-2.5)),
        axis.line = element_line(size=axisLineSize, colour=axisLineColour))
detect.respondToIdiot.nomoral_expt5

ggsave("img/modelRespondToIdiot_detect_nomoral_expt5.png", detect.respondToIdiot.nomoral_expt5, height=8, width=7.5)  


detectplots_nomoral_expt5 <- plot_grid(detect.respondToIdiot.nomoral_expt5 + ggtitle("Random Opponent") + theme(plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold")),
                         detect.recurse.nomoral_expt5 + ggtitle("Recursive") + theme(plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold")), 
                         detect.results_expt4 + ggtitle("Red Utility") + theme(plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold")),
                         detect.results_expt5 + ggtitle("Blue Utility") + theme(plot.margin = marginDim,
        plot.title = element_text(size=plotTitleSize, hjust = 0.5, face = "bold")), 
                         nrow=1,
                         rel_widths=c(10.72, 10, 10, 10))


detectrow_nomoral_expt5 <- plot_grid(detectplots_nomoral_expt5, detectorLabel, nrow=1, rel_widths=c(96, 4))

withLegend_nomoral_expt5 <- ggdraw() +
  draw_plot(detectrow_nomoral_expt5) +
  draw_plot(legend, x=0.415, y=0.32) 
ggsave("img/detectPredictions_expt5.png", 
       withLegend_nomoral_expt5, 
       units="cm", height=4.7, width=17.8)